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ABSTRACT 


This  thesis  is  concerned  with  numerical  ap¬ 
proximations  to  the  solutions  of  differential  equations. 

The  elements  of  the  finite-difference  calculus 
are  presented,  and  employed  to  develop  several  useful 
formulas.  Elementary  methods  for  the  solution  of  differ¬ 
ential  equations  are  presented  first,  followed  Toy  more 
sophisticated  techniques.  Instability  of  difference 
approximations  to  initial  value  problems  is  discussed  at 
some  length,  and  an  algorithm  which  attempts  to  minimize 
this  difficulty  is  presented. 

The  classification  of  partial  differential 
equations  into  hyperbolic,  parabolic,  and  elliptic  types 
is  shown,  and  methods  for  the  solution  of  all  three  types 
are  given. 


A  discussion  is  given,  approaching  the  Green's 
function  from  a  numerical  point  of  view.  It  is  shown  that 
this  perspective  leads  to  new  insight  into  the  subject.  An 
analytic  expression  for  the  Green's  function  of  a  large 
class  of  differential  operators  is  derived  by  an  essential¬ 
ly  numerical  technique. 
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CHAPTER  I 
Introduction 

Section  1.1  A  survey  of  the  Thesis 

One  of  the  important  problems  common  to  all  of 
the  physical  sciences  is  the  formulation  and  solution  of 
differential  equations.  A  large,  and  still  rapidly  growing, 
literature  has  been  developed  in  reponse  to  this  need. 

As  might  be  expected,  the  passage  of  time  has 
seen  many  of  the  less  difficult  problems  solved,  and  so, 
the  problems  now  posed  are  often  quite  difficult.  More  and 
more  often,  new  problems  are  intractable  by  traditional 
techniques.  For  instance,  it  is  known  that  no  expression 
in  terms  of  elementary  functions  can  be  given  for  the 
solutions  of  the  equations 

y'  -  x2  +  y2,  (l.l.l) 

and 

y"  =  6y2  +  x. 

Furthermore,  even  if  an  expression  Is  available  for  the 
solution  of  a  differential  equation,  its  evaluation  may 
prove  quite  tedious;  an  example  of  this  situation  Is  given 
in  Chapter  III. 

For  many  years,  numerical  methods  for  solving 


•  <  '■  ■  W.  ■  ;  '."’1  :  ,  3.  I- 


o  :Vr-  ■  ;  ,jJ:  n-<  \3  3.x...  ul  '■  /  :  •  t:33:vxx;9xrvX£>-'3 

;  3.  O  .  3  •  'VO-  ' 


differential  systems  In  an  approximate  sense  have  been 
known,  but,  on  account  of  the  large  amounts  of  arithmetic 
involved  in  carrying  them  out,  they  have  not  been  extensively 
employed  until  recently. 

Since  the  early  and  middle  1940' s,  electronic 
digital  computers  have  become  available.  Currently  available 
machinery  can  perform  upwards  of  250,000  arithmetic  operations 
per  second,  and  at  a  relatively  modest  cost  (something  like 
.01/  per  operation  for  an  IBM  7090) .  Thus,  methods  of 
numerical  approximation  have  been  rendered  practical,  and 
therefore  a  great  deal  of  recent  research  has  been  directed 
toward  the  further  development  of  these  techniques. 

This  thesis  attempts  to  review  the  major  techniques 
developed  for  the  numerical  integration  of  differential,  and 
especially  ordinary  differential,  equations.  It  does  not 
pretend  to  be  complete,  since,  to  do  so,  it  would  have  to 
be  many  times  its  present  size;  neither  does  it  go  into  the 
detailed  development  of  a  topic,  unless  this  is  in  some  way 
revealing  of  general  techniques. 

The  second  chapter  presents  the  basic  results  of 
the  finite-difference  calculus,  and  compares  it  with  the 
more  familiar  differential  calculus.  Several  results  are 


derived  for  later  use  in  the  thesis. 
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Chapter  III  attempts,  by  way  of  elementary 
techniques  and  examples  to  provide  a  taste  of  the  numerical 
approximation  of  differential  equations.  The  purpose  is  to 
exhibit  ideas,  free  from  too  much  obscuring  detail.  Elemen¬ 
tary  approximations  to  initial  and  boundary  value  ordinary 
differential  equations,  and  to  hyperbolic,  parabolic,  and 
elliptic  partial  differential  equations  are  presented  and 
discussed.  The  problem  of  stability  is  illustrated  by  means 
of  numerical  example. 

In  Chapter  IV,  practical,  and  more  advanced 
methods  are  discussed,  for  ordinary  differential  systems  only. 
The  important  multi-step  integration  techniques  are 
developed,  and  their  stability  properties  are  discussed  in 
some  detail.  Runge-Kutta  methods  are  also  discussed,  and 
programs  for  multi-step  and  Runge-Kutta  integrations  are 
given.  Methods  for  two-point  boundary-value  problems  are 
mentioned  briefly. 

In  Chapter  V,  we  pursue  a  line  of  thought  partially 
developed  in  Chapters  II  and  III.  We  consider  the  matrix 
equation  obtained  by  replacing  a  differential  equation  by 
a  finite-difference  approximation.  An  expression  is  obtained 
for  the  solution  of  the  algebraic  problem,  and  it  is  shown 
that,  in  the  limit  as  the  tabular  interval  approaches  zero. 
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the  finite-difference  approximation  yields  an  exact 
solution  in  terms  of  Green's  functions.  Further  an  expression 
for  the  Green's  function  of  a  quite  general  differential 
operator  is  obtained  by  this  method.  It  is  believed  that 
some  of  the  results  of  this  chapter  are  new. 

The  final  chapter  summarizes  the  major  results 
of  the  thesis,  and  indicates  some  topics  on  which  current 
research  is  proceeding. 

Each  chapter  is  preceded  by  a  brief  synopsis. 

Section  1.2  Notation  Used  in  the  Thesis 

The  notation  used  is  largely  standard,  the  major 
deviation  occurring  in  the  programs  of  Chapters  III  and  IV. 
These  programs  have  been  written  using  a  subset  of  an 
algorithmic  language  proposed  by  Iverson  [27]* 

According  to  this  notation,  ordinary  variables 
are  written  as  lower-case;  vectors  as  underscored,  lower¬ 
case;  and  matrices  as  underscored,  upper-case  Roman  characters. 
An  element  of  a  vector  is  indicated  by  a  subscript,  without 
the  underscore;  a  column  vector  of  a  matrix  is  indicated 
by  a  subscript,  and  the  underscore  is  retained. 

Each  step  of  a  program  is  numbered,  for  possible 
further  reference.  The  notation 
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implies  that  the  variable  a  receives  the  value  of  the 
expression  (3.  The  form 

u.pj  (r^r^,  .  .  . )  ►( s^jSgj  .  .  . ) 

means  that  a  and  |3  are  compared,  and  that  if  a  stands  in 
the  relation  r^  to  p,  control  passes  to  step  s1;  if  in 
relation  r^,  to  s^,  etc.  The  relations  are  of  the  form 
<_,  >j  />  and  so  on. 
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CHAPTER  II 


THE  CALCULUS  OF  FINITE  DIFFERENCES 

This  chapter  introduces  the  finite  difference 
calculus.  It  attempts  to  put  it  in  a  perspective  which 
makes  it  seem  as  natural  a  descriptive  device  as  differential 
calculus . 


The  elementary  finite  difference  shifting, 
differencing,  and  averaging  operators  are  defined,  and  some 
formal  identities  are  obtained. 

Expansions  of  the  familiar  differentiation 
operator  in  terms  of  differencing  operators  are  derived. 

Some  formal  analogies  between  finite  difference 
and  differential  calculi  are  shown.  In  particular,  an 
example  is  given  of  the  solution  of  a  difference  equation  by 
techniques  borrowed  directly  from  differential  calculus. 


Section  2.1  Introduction 

In  the  elementary  study  of  differential  calculus, 
we  may  approach  the  concept  of  a  derivative  in  the  following 
way : 

We  have  an  intuitive  notion  of  what  is  meant  by 
the  "slope"  of  a  given  curve  f(x)  at  any  given  point  x  .  We 
can  easily  speak  of  the  "average  slope,"  S(x)  of  f(x)  in  the 
interval,  say,  xq<x  <fx^,  defined  as 


X  +xn  f  (x-,  )  -  f(x  ) 

S(V>  -  ^ 

d  X1  '  xo 


(2.1.1) 


By,  "the  slope  of  f(x)  at  xo"  we  then  understand 
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the  limit  (it  it  exists)  of  (2.1.1)  as  x,  is  made  to 
approach  xQ.  This  limit  is  termed  the  derivative  of  f(x) 
at  x  =  x  .  The  derivative  of  f(x)  at  any  x  is  written 


f 5  (x) 


df  (x) 

dx 


Df(x) 


lim 
X-j-*  x 


f(xx)  -  f ( x ) 

X^-  X 


(2.1.2) 


It  is  worth  noting,  in  passing,  that  this 
definition  requires  at  least  the  continuity  of  f(x)  at  x, 
and  that,  in  the  physical  world,  this  condition  is  seldom 
fulfilled,  due  to  the  atomicity  of  most  material  objects. 
Thus,  in  the  real  world,  this  concept  is  usually  an 
idealisation,  since  the  limit  of  (2.1.1)  cannot  actually 
be  taken. 


As  an  example  of  this,  the  derivation  of  the 
well  known  string  equation  proceeds  by  considering  the 
effect  of  elastic  forces  on  a  very  short  section  of  string, 
applying  Newton’s  second  law  of  motion,  and  deriving  an 
equation.  The  length  of  the  section  of  string  is  then 
reduced  without  limit,  and  a  differential  equation  results. 
Really,  this  process  is  not  valid,  since  strings  are  known 
to  be  composed  of  atoms  and  molecules  of  non-zero  size. 
However,  for  practical  purposes,  i.e.,  for  the  sake  of  a 
manageable  physical  model,  this  distinction  is  trivial. 
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Nevertheless,  we  have  illustrated  here  the 
important  fact  that  most  of  the  differential  equations 
used  in  the  physical  sciences  are  derived  by  obtaining  an 
equation  at  separated  points  and  then  letting  the  points 
approach  one  another.  This  thesis  will  explore  some  of 
the  implications  of  this  procedure,  especially  in  Chapter  V 

Again  in  passing,  it  is  to  be  admitted  that 
some  differential  equations  are  not  obtained  in  this  way. 
For  instance,  in  quantum  mechanics,  Schr<Sdinger 1  s  wave 
equation  is  obtained  essentially  by  inspection,  experiment, 
and  inspiration.  This  is  also  the  case  for  Charles1  law 
for  gases  and  Newton’s  laws  of  motion. 

Let  us  now  go  on  to  the  finite -difference 
calculus  itself. 

Section  2.2  Finite  Difference  Operators 

We  consider  a  function,  say,  f(x),  defined  at 
the  abscissae 

x  -  xq,x1,x2, . . . ,xn, 

where,  for  some  real  number  h  (usually,  but  not  necessarily 
thought  of  as  being  small), 

x .  =  x  +  ih.  :  i=0,l,2,.»*.,n. 
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We  will  usually  write 

f (x. )  =  f±. 

A  number  of  operations  can  now  be  defined  on  the 
f. .  The  names  of  the  operations,  the  symbols  used  to 
represent  them,  and  the  effect  of  each  operator,  acting  on 
f,  follow: 

Forward  Shift  operator,  E  ,  Ef^  =  fi+1  ^  2  l) 

Forward  Difference  operator,  A  ,  Af^=f^+1~f^  (2.2.2) 

Backward  Difference  operator,  V  ,  Vf^=  f^-f^^  (2.2.3) 

Central  Difference  operator,  5  ,  Sf^  =  ^  j_+]_/2~^i-l/2 

(2.2.4) 

Averaging  operator,  p.  ,  \if±  =  l/2(f i+]_/2+f  i-l/2^ 

(2.2.5) 

Differentiation  operator,  D  ,  Df.  =  (d^^ ) 

1  ax  x-x.  (2.2.6) 

Superscripts  are  used  to  indicate  repeated 
application  of  an  operator.  That  is,  representing  any  of 
the  above  operators  by  0, 

0nf±  -  o(on_1fi) .  (2.2.7) 

For  example. 


Z2fi  =  E(Ef.)  =  E(f.+1)  =  f1+2, 


(2.2.8) 
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■  .  ■  ' 


.  ;  i,  ,  ..  .  ■■  V  .  ,  w-  A 


- 

,  '  "'(■  .  ■  ■  r/  ;qo  s  ■■■■■$£  ©r l 1 
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and 


f 


i+n 


(2.2.9) 


We  will  treat  these  operators  formally,  to  obtain 
some  useful  identities.  The  rigorous  justification  for  most 
of  these  identities  can  only  be  given  a  posteriori,  and  in 
the  context  of  a  particular  function,  f(x),  or  a  particular 
class  of  functions. 


First,  we  will  express  the  shift  operator,  E,  in 
terms  of  the  differentiation  operator,  D.  We  have 


Ef.  = 

i 


fi+lj 


that  is, 

Ef(x1)  =  f(x1+h).  (2.2.10) 


But,  by  Taylor* s  theorem, 

h2D2f (x. ) 

f(xi+h.)  =  f(x.)  +  h.Df(xi)  +  — 2! - —  +  •••  •  (2.2.11) 


That  is, 

f  ( x .  +h. ) 
v  1  ' 


( 1+h.D- 


h2D2 

2! 


+ 


+...)f(xx) 


(2.2.12) 


We  recognize  the  series  in  (2.2.12)  as  the  formal  series  for 
the  exponential  function.  Hence,  we  write 


Ef.  =  f  ( x .  +h. )  =  ehDf(x.) 
1  x  1  '  v 


(2.2.13) 


o::  “  1  X ;  •  .  i  j  .  a^tol  v  j'-;  aamvj  Lit 

ctoifsoUllSaut  &-,otio'(gicx  axtfE  JStfrrsibi:  Jjj'iasu  sirros 

.  i.  i  '  ::  [  r  :a  9..  .  i:  'BO  ■.  jj::i ...  -  :  >a  :  •  ../  lo 

. -  • 

'  r  •'  ,  ■!  )'l  ,  a-:  iX  nu  .  n&lwjj  :  '  3  39  vJn'vO  aria 


■  c  I va  .  a  ^ 


••  (  .  '  ■  V" 


. 


‘  .  2 ) 


.  38  iJE.  r  I  i  8.8  (21,2  2')  ft./  :SJ.iigOO£K).  sW 


■  *j  03J./1  r  j*r>  Jf.... 9  3ii.t 

t  r  )l  -  rl3 
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and  thus,  we  obtain  the  formal  relation 


T7  h.D 

E  =  e 

(2.2.14) 

Now,  we  will  express  the  three  difference 

operators 

and  the  averaging  operator  in  terms  of  E. 


We  have 

Af.  =  f.  -  f.  =  Ef.  -  f.  =  ( E-l ) f . , 

l  l+l  i  l  l  k  l 

(2.2.15) 

and  hence,  formally. 


A  =  E-l. 

(2.2.16) 

Now, 

Vf.  -  f.-f.  ,  =  f . -E-1f .  =  ( 1 -E_1 ) f . , 

l  l  l-l  l  i  v  '  l 

(2.2.17) 

and,  therefore. 

V  =  (l-E-1 ) . 

(2.2.18) 

We  note  that  E  1  is  a  "backward  shift  operator." 


Also, 

6fl  =  fi+l/2  -  fl-l/2  =  El/2fi  -  E'l/2fl 

-  (E1/2-E“1/2)fi, 

(2.2.19) 

and  thus, 

5  =  (e1//2-E-1/2)  . 

(2.2.20) 

:  .  '  .  .  ) 


'  ' 


ri  =  .o  .■  '  ■'  ■ "  •  .> :•$)  -  iJttf  ■  |  j 


..  ;j  ’  ■  ,  ;o. ..  '  .r;  ..  r&  oi,.  bn...-: 


i.'A 


Y-CXi  .  eon • 


h 


■.  - ;  -  •  ■  i-'j.  y  i"  2  J;  7.  rsrfj-  o a  9W 

,  .  ;  a 


i1""-'  -  .  X  3  ■  a,..-,'  -  ' :\,  Y"  3 


I 


,  ':0a  :  ,  ■  .  . 
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Finally, 


^  =  l/2(f1+l/2+fi.l/2)  =  |(E1/2+E-V2)fi, 

(2.2.21) 

and  so 

n  =  |(E1/2+e-1/2). 

(2.2.22) 

We  can  combine  these  relations  with  (2.2.14)  to 
obtain  expressions  in  terms  of  D. 

We  have  : 


A  =  E-l  =  ehD-l ; 

(2.2.23) 

V  =  1-E-1  =  1  -e_hDj 

5  =  (E^+E-l/2)  =  (ehD/2_e-h.D/2):=2  sinh.  h|. 

(2.2.24) 

(2.2.25) 

and 


H  =  |(E1/2+e-1/2)  ,  l(ehD/2+e-hD/2)  =  cosh  hD  _ 

(2.2.26) 

A  further  relation  which,  we  will  need  later  can 
now  be  derived.  Substituting  (2.2.25)  and  (2.2.26)  into 
the  identity 


cosh.^  ^-h.D  -  sinh.^  ^h.D  =  1, 

(2.2.27) 

we  obtain. 

2  1  =2  , 

\i  -  If  5  =  1* 

(2.2.28) 

1'e"  m.  =  (1  +  }  52)1/2. 

(2.2.29) 

!  '  '•  - ' 


.  ■-  V-  . 


. 


.  .  !  '  J 


v 


. 

'  ■’  ■ ■  ;  1'; .  ,  ■ :  £  ;  ■.  ■  '  .  :  v  •'  ■xu  -J  roa 


,  f  -  '  i  M  X  L ..  -  Cf  Cf:  "  ff  30'- 
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In  fact,  each  of  the  operators,  (2.2.1)  to  (2.2.6), 
can  be  expressed  in  terms  of  each  of  the  others  [l],  but 
the  relations  derived  above  are  representative,  and  will 
suffice  for  our  purposes. 

Section  2.3  Derivatives  in  Terms  of  Difference  Operators 

From  the  point  of  view  of  solving  differential 
equations,  it  is  very  useful  to  be  able  to  express  the 
operator  D,  and  its  powers,  as  evaluable  expressions  of 
difference  operators.  Our  motivation  in  this  section  is 
partly  formal  and  theoretical,  and  partly  practical. 

Writing  (2.2.23)  as 

ehD  =  1  +  A,  (2.3.1) 

and  taking  natural  logarithms  of  both  sides,  we  obtain 

h.D  =  1  n  ( 1+A )  .  (2.3.2) 

This  is,  itself,  not  an  evaluable  expression  in 
terms  of  A,  but  we  can  make  it  so  by  expressing  the  right-hand 
side  as  a  Taylor  series.  Then, 

hD  =  In (1+A)  =  (A-|A2+iA3+. . . ) .  (2.3*3) 

In  fact,  for  derivatives  of  any  order,  m,  we  have 


, 

S9;::CK'i*iLrq  -1.,/D  ol  ioxllj,];:.! 


'  ■:•■■■  :.,i  gt  :  77. .gig  y  : 7.  ?ffc- .tJ-.-d/jps 


: 


: 


7.'.  7'  G  •  .7  0  :  ,  l  '  977  7  ,  .  •  G.GiiT  - 

.  1  :i3  <  9  -■<  ■  O-  ,  .1  .  *•-,  f:  iC  ■  io  ■">  (T/‘ J  -• 


■ 


.  -G  ■':  (  '  7 '  ■'(  ::  W  .  ■  .  .7 


(2.3.4) 


hmDm  =  (A  -  t  ^  A-5 


In  particular,  for  m  =  2,  [2] , 


b2D2  =  -  AJ  +  A 


3  j.  li  a4  .  |a5  +...). 


(2.3.5) 


In  a  similar  fashion, 


hmDm  =  (V  +  |  V2  +  j  V3  +. . . ) 


m 


(2.3.6) 


Now.  rewrite  (2.2.25)  as. 


hD  =  2  sinb  2  5/2 


(2.3.7) 


i .  e  .  , 


/h.D\m  /sinh^S/^N 
'  5  “  '  572  ~ ^ 


m 


(2.3.8) 


-l 


The  Taylor  series  for  sinh.  x  is,  [3], 


-1  1  3  ^  S  h  7 

sinh.  x  =  x  -  -g-  x^  +  ^  x^  -  x  + 


(2.3.9) 


Hence , 

n  m^m  Km/n  6  2  ,  3  c4  5  c6  \lrl 

hD  =  5  (1“  2^  +  OTO  6  “  14733F  5  . . . ) 


(2.3.10) 


■ 


' 
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In  particular,  for  m  =  2, 


2  4 

,  2n2  r2/,  5  ,6^ 

h  D  =  5  (1  -  Y2  +  90  - 


560 


+•••)• 


(2.3.11) 


We  notice  that  the  coefficients  in  (2.3.10)  become 
small  more  rapidly  than  those  of  (2.3.4)  or  (2.3.6).  That 
is,  (2.3.10)  tends  to  converge  more  rapidly  than  the  others. 
This  seems  to  be  a  typical  feature  of  central  difference 
formulae  in  general. 

One  objectionable  feature  of  ( 2.3.10)  is  that  for 

odd  values  of  m,  the  right-hand  side  contains  only  odd 

powers  of  5.  Hence,  in  order  to  use  this  formula,  we  would 

1  3 

need  to  know  our  function  f(x)  at  the  points  xq+  75-h.,  XQ  +  2 
and  so  on.  To  see  this,  we  write  (2.2.20): 

5  =  E1/2  -  E-1/2  =  E_1/2(E-l) .  (2.3.12) 

Therefore , 

5n  =  E~n/2(E  -  l)n. 

Clearly,  for  n  odd,  this  series  contains  a  half 
power  of  E  in  every  term,  and  implies  the  need  for  values  of 
f  at  the  mid -tabular  points. 


To  circumvent  this  difficulty,  we  multiply  the 


■  . 


'  •  ■  :•  .  . 


■«.  ;  sa';. 


■  '  WQ 


■  a  .  .. 
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right-hand  side  of  (2.3.10),  by  to  obtain 


,  m^m  mKm 
h,  D  =  |jl  5 


(a  1(l  - 


6  +  3 

2IT  OTO 


A 

o 


m 


(2.3.13) 


According  to  (2.2.29). 


-1 


(1  +  £ 


i  5  3  .4 

1  3“  12o  °  •  •  • 


(2. 3.14) 


Thus,  multiplying  the  two  series  together,  we 


obtain. 


h.  D 


=  ([l6)m(l  -  |  62  +  U  54-...)"1 


(2.3.15 


This  series  involves  only  the  tabular  points  of 
f(x),  for  even  or  odd  m.  It  could  be  used  in  place  of 
(2.3.10),  even  for  even  m,  except  for  its  somewhat  inferior 
convergence . 

Relations  (2.3.10)  and  (2.3.15)*  as  well  as  other 
generally  useful  derivative  formulae,  are  tabulated  in  [2], 
for  various  values  of  m. 


Section  2.4  Comparisons  of  Finite  Difference  with 

Differential  Calculus 


We  can  again  consider  the  "average  slope" 


' 


‘ :  i  i 


.  n,  B|cfn 


stfnioq  2 ;  erJ-:! 


. 
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formula  (2.1.1): 


x  +xn  f  (x,  )  -  f  (x  ) 

S(^l  )  -  1 


X-,  -X 

1  o 


More  generally. 


x.+x.,n  f(x.,n)-  f(x.) 

st  1  y+1)  =  1+1 — 

2  x.+1-  x. 


Af . 
1 

Ax. 

l 


Thus , 


f  ’  (x) 

v  'X=X. 

1 


n  .  Af  . 

nf.  lim  1 

i  Ax.— ►  0  Ax.  ’ 

l  l 


and  so,  for  small  Ax1  =  h,  we  expect  that 


Df .  %  Afi 
1  h. 


(2.4.1) 


(2.4.2) 


or,  in  other  words,  that 

h.D  %  A . 


(2.4.3) 


We  are  thus  led  to  suspect  that  the  differentiation 
operator,  D,  and  the  various  difference  operators,  might, 
formally  at  least,  behave  similarly.  In  fact,  they  do  so 
as  we  shall  show,  although  there  are  differences  which 


IJ  ■•.o'l  ,oa  bnis 


. 
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prevent  the  analogy  from  being  complete. 

For  instance,  for  suitable  functions  f  and  g,  and 
constants,  A  and  B, 

D(Af  +  Bg)  =  ADf  +  BDg, 

and 

A(Af1  +  Bg;.)  =  AAf±+  BAg± .  (2.4.4) 


In  this  case,  the  analogy  is  complete.  In  the 
following,  we  consider  derivatives  and  differences  of 
quotients  and  products,  and  the  analogy  is  somewhat  damaged. 

In  differential  calculus, 

D(f-g)  =  f.(Dg)  +  g-(Df).  (2.4.5) 

In  difference  calculus, 

A(f.g)  =  (E-l)(f-g)  =  E(f.g)  -  f-g.  (2.4.6) 

Now,  by  definition, 


E(f-g)  =  f(x+b)g(x+h)  =  (Ef)-(Eg). 


(2.4.7) 


, 


(  ■>:  IV 
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Hence , 

A(f.g)  =  (Ef).(Eg)  -  f.g 

=  (Ef)-(Eg)  +  (Ef)-g  -  (Ef)-g  +  f-g 
=  (Ef)-(Eg-g)  +  (Ef-f)-g 
=  (Ef) . (Ag)  +  g. (Af) .  (2.4.8) 

By  a  similar  argument, 


A(f-g)  =  f  • (Ag)  +  (Eg) • (Af )  .  (2.4.9) 

Taking  the  average  of  (2.4.8)  and  (2.4.9)?  we  can 
obtain  the  symmetric  form: 

A(f-g)  =  f *  ( Ag)  +  (^g.(Af)  (2.4.10) 

The  corresponding  expressions  for  the  backward 
difference  operator  are : 


V(f.g)  =  (E-1f ) * (Vg)  +  g-(Vf)  (2.4.11) 

=  f-(Vg)  +  (E-1g) * (Vf )  (2.4.12) 

=  (E_bl)  f.(Vg)  +  gU)  g.(Vf). 

(2.4.13) 


20  - 


All  of  equations  (2.4.8)  to  (2.4.13)  are  formally 
similar  to  (2.4.5)  except  for  the  operators  involving  E. 

We  also  note  that  for  h.  — ►  0,  all  of  these  operators  involving 
E  tend  to  unity. 

Again  in  differential  calculus, 

D(f/g)  =  K-(Pf)-  f-(PR) 

g2  (2.4.14) 

whereas,  in  difference  calculus, 

A (f/g)  =  (E-l )  ( f/g)  =  E(f/g)  -  f/g 

=  (Ef)/(Eg)  -  f/g  =  g^Ef')  ~ 

g* (Eg) 

=  g • (Ef )  -  g • f  +  g-f  ~  f * (Eg) 

-  g • (Eg) 

=  g-(Ef-f)  -  (Eg-g).f 

g* (Eg) 

=  g-(Af)  -  f.(Ag)  .  (2.4.15) 

g  * (Eg) 

g • ( Vf )  -  f  * ( Vg) 

g*E_1g 


Similarly, 

V(f/g) 


(2.4.16) 
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Once  again,  (2.4.14)  is  formally  the  same  as 
(2.4.15)  and  (2.4.16)  except  for  the  operators  E  and  E 

In  differential  calculus,  the  operation  of 
integration  is  inverse  to  D.  In  difference  calculus, 
summation  is  inverse  to  differencing,  as  we  shall  now  show. 

Consider  the  definite  summation 
b 

V  g(x)  =  g(a)  +  g(a+h)  +  g(a+2h)+. . .+g(b) ,  (2.4.17) 

a 

where  g  is  a  function,  and  a  and  b  are  real  numbers  such,  that 
a<b,  and  (b-a)  =  nh. 

Suppose  that  we  know  some  function  such  that 

Af(x)  =  f(x+h)  -  f(x)  =  g(x) .  (2.4.18) 

In  other  words,  f  is  an  "anti-difference"  of  g, 
which  we  write 

f(x)  =  &_1g(x) .  (2.4.19) 


Substituting  (2.4.18)  into  (2.4.17)*  we  have. 


’ 


. 
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b 


g(x)  =  g(a)  +  g(a+h)  +  g(a+2h)+...+  g(b-h)+  g(b) 


a 


f(a+h)  -  f(a) 


+ 


f(a+2b)  -  f(a+b) 


+ 


f (a+3h) -f (a+2h) 


+. . .+  f  (b)  -  f (b-h) 


+  f (b+h) -f (b ) 


(2.4.20) 


We  observe  that  all  the  terms  on  the  right  in 
(2.4.20)  cancel,  except  for  -f(a)  and  f (b+h) .  Hence, 


b 

=  f(b+h)  -  f(a). 
a 


(2.4.21) 


This  corresponds  closely  to  the  "fundamental 
theorem"  of  the  integral  calculus : 


x)  dx 


f(b) 


f(a),  where  Df(x)  =  g(x) 


(2.4.22) 


In  fact,  the  method  for  obtaining  (2.4.21)  suggests 
a  method  of  deriving  (2.4.22).  For,  consider  the  definition 
of  the  definite  integral, 

b  .  b  .  b 

J  g(x)dx  =hi^Q  ^hg(x)  =h^1Q  h^g(x)  (2.4.23) 

a'  a 


1 r  \  ,  - .. 1 : ) 

■ 


■  i;;  "o,  f‘rae‘rv*3cUf 


-  - 


-• 


-  23 


Also, 

Df(x)  =  llm  Af (x)  . 
h— ►  0  h 

Furthermore,  as  h.-*0, 

f  (b+h)  f  (b) . 

The  above  argument  does  not,  of  course,  even  begin 
to  be  rigorous.  It  does,  however  provide  some  additional 
motivation  for  the  "fundamental  theorem"  beyond  that  usually 
given  in  elementary  calculus  courses [6] j  that  is,  in  a  sense, 
it  shows  "why"  the  theorem  is  true. 

There  are  many  more  formal  similarities  between 
the  difference  and  differential  calculi  [4  ]  .  We  consider 
briefly  difference  and  differential  equations  in  the 
following  section. 

Section  2.5  Finite  Difference  Equations 

In  the  differential  calculus,  the  general  form  of 
an  ordinary  differential  equation  of  order  n  is 

F(Dnf,Dn_1f, . . . ,Df,f,x)  =  0  (2.5.1) 

where  F  is  any  function  of  n  +  2  arguments,  and  f(x)  is 


. 


■ 
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an  unknown  function. 

Corresponding  to  (2.5*l) >  the  general  difference 
equation  of  order  n  may  he  written 


F(Anf,An  1f , . . . ,Af,f ,x)  =  0. 


(2.5.2) 


Since  we  can  always  write  difference  operations 
in  terms  of  shift  operations,  as 

A^vf  =  a^Ekf  +  a^^E^  1f+.  •  . +a^Ef  +  aQf ,  (2.5*3) 


it  is  therefore  clear  that  an  alternate  expression  for 
(2.5*2)  is 

G(Enf,En-1f, . . . ,Ef,f,x)  =  0,  (2.5*4) 


or  also. 


G  f  ( x+nh. ) , 


f (x+n-lh) , . . . ,  f(x+h),  f(x),x 


0.  (2.5.5) 


In  the  sequel,  purely  for  simplicity  of  notation, 
we  will  assume  h  =  1.  It  is  clear  that  appropriate 
transformations  can  he  made  to  difference  equations  with 
arbitrary  h.,  to  bring  them  to  this  form. 
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As  we  shall  presently  demonstrate,  many  of  the 
elementary  tricks  and  techniques  for  solving  differential 
equations  can  be  applied,  with  slight  modifications,  to 
difference  equations. 

For  example,  consider  the  difference  equation 

f ( x+2 )  -  7  f(x+l)  +  10f(x)  =  0.  (2.5.5) 

As  a  trial  solution  of  (2.5.5.)?  we  try  f(x)  =  ax. 
Substituting  this  trial  function  into  the  equation,  we  have 

ax+2  -  7ax+1  +  10ax  =0.  (2.5.6) 

i.e.,  ax(a2  -  7a  +  10)  =  0  (2.5.6) 

Now,  the  case  of  a  =  0  is  a  trivial  solution,  and 
we  disregard  it.  The  remainder  of  (2.5.6)  is  a  quadratic  in 
a,  and  has  roots 

a  =  2  and  a  =  5 • 


is 


Hence,  it  seems  that  a  general  solution  of  (2.5*5) 


f(x)  =  A-2X  +  B-5X. 


(2.5.7) 


. 
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That  (2.5.7)  is  indeed  a  solution  may  be  verified 
by  direct  substitution  into  (2.5.5). 

However,  it  can  also  easily  be  verified  that 

f (x)  =  Au(x)2X  +  Bv(x)5X  +  w(x)  (2.5.8) 

is  also  a  solution,  for  any  u(x)  that  has  some  constant 
value,  for  any  v(x)  that  has  a  constant  value,  and  for  any 
w(x)  that  is  zero,  at  all  the  tabular  points. 

For  instance, 

f(x)  =  A2Xcos(2ttx)  +  b5Xtan(-Trx  +  ^)  +  sinvx  (2.5.9) 

is  also  a  solution  of  (2.5*5)*  since  cos(2t rx)  is  1, 
tan(7rx+^-)  is  ,  and  sin-rrx  is  zero  at  all  o  the  tabular 

points,  x  =  0,1,2,... 

If,  then,  we  are  concerned  only  with,  the  solution 
at  tabular  points,  as  is  indeed  the  case  in  this  thesis,  then 
(2.5.7)  provides  an  unique  solution  when  A  and  B  have  been 
fixed . 

Clearly,  in  order  to  determine  the  constants 
A  and  B,  we  require  two  additional  conditions.  This 
requirement  is  entirely  analogous  to  the  two  conditions 
needed  for  a  linear,  second  order,  differential  equation. 


' 


■ 
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These  conditions  may  take  the  form  of  one 
equation  at  each  of  two  points  (boundary  conditions),  or 
two  independent  equations  at  one  point  (initial  conditions). 

For  example,  we  might  have  boundary  conditions  like 
f(l)  =  12,  and  f(3)  =  258.  (2.5- 10) 


Hence,  we  could  write  the  linear  equations: 

12  =  2A  +  5B 

258  =  8A  +  125B,  (2.5.H) 

which,  we  can  readily  solve,  to  find, 

A  =  1,  B  =  2, 

and  hence. 


f (x)  =  2X  +  2- (5*) 


x- 


(2.5.12) 


like , 


have 


Alternatively,  we  might  have  had  initial  conditions, 
f(l)  =  9,  and  Af ( 1 )  -  24.  (2.5-13) 

Substituting  the  first  of  these  into  (2.5*7)?  we 

9  =  2A  +  5B.  (2.5*14) 


.  i  •.  .'0  i  'Jj  G  .//-■  J  .  ■;  >  '  rt^upe 
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Differencing  (2.5.7) >  we  have 
Af (x)  =  AA ( 2X )  +  B  A(5X) 

=  A(2X+1-2X)  +  B(5X+1-5X) 

=  2XA  +  4B5X.  (2.5.15) 

Substituting  the  second  of  (2.5-13)  into  (2.5*15) * 

we  obtain 

24  =  2A  +  20B .  (2.5.16) 


Prom  (2.5.1^)  and  ( 2.5.16),  we  have,  immediately, 

A  =  2,  and  B  =  1, 

and  hence, 

f(x)  =  2 • 2X  +  5X  =  2X+1  +  5X-  (2.5.17) 

More  difficult  problems  require  additional  results 
about  antidifferences,  and,  in  particular,  about  factorial 
polynomials  and  Stirling  numbers  [5]*  However,  the  preceeding 
example  will  serve  to  illustrate  the  basic  approach. 
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CHAPTER  III 


ELEMENTARY  METHODS  FOR  THE  NUMERICAL  SOLUTION  OF 
DIFFERENTIAL  EQUATIONS 


This  chapter  employs  the  techniques  of  the  previous 
chapter  -  and  some  others  -  for  the  solution  of  differential 
equations . 


The  Taylor  series  and  Euler  approximation  methods 
for  sets  of  first  order  initial  value  differential  equations 
are  developed.  The  important  Cauchy  existence  and  uniqueness 
theorem  rising  out  of  the  Euler  approximation  is  stated  and 
discussed,  and  stability  of  initial  value  problems  is 
illustrated  by  means  of  an  example  of  an  unstable  approximation. 

Boundary  value  ordinary  differential  equations, 
including  eigenvalue  problems,  are  briefly  introduced. 

The  general,  second-order,  quasi-linear  partial 
differential  equation  is  classified  into  hyperbolic,  parabolic, 
and  elliptic  types.  The  concept  of characteristic  curves 
is  introduced  and  the  solution  of  hyperbolic  equations  by 
the  method  of  characteristics  is  described  in  detail. 

Simplified  outlines  are  given  of  methods  for  solving  parabolic 
and  elliptic  systems. 


Section  3-1  Introduction 

Increasingly,  the  new  differential  equations 
which  arise  in  practice,  do  not  admit  of  a  solution  in  closed 
form.  Even  if  they  do,  the  labour  involved  in  tabulating 
the  solution  may  be  quite  large. 

An  example  given  in  Modern  Computing  Methods  [11] 


is  illustrative: 
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The  solution  of  the  simple  equation 


dy 

dx 


2 3L 


1  — x 


4 


0, 


(3.1.1) 


can  be  shown  to  be 


y 


A  t  l-  "1  ' 

=  Mirp  exp  ( tan  x, 


(3.1.2) 


The  very  slightly  more  complicated  equation 


dy  _  2y 


dx 


1  — x 


4 


=  x 


(3.1.3) 


has  the  solution 


/1+X\  ^ ^  -1  ^  P  /1-X\  ' 

y  =  (^Tp  exp  ( tan  x)  /x (tzt)  exP 


1/2 
-X'  / 

1+X‘ 


(-  tan  1x)  dx  +  A 


(3.1.4) 


To  evaluate  (3.1.4)  is  no  mean  task;  in  particular, 
exponential,  arctangent,  and  square  root  functions  must  be 
evaluated,  and  a  numerical  quadrature  performed. 

The  present  chapter  presents  an  alternative 
procedure  -  direct  numerical  approximation  from  the  differential 
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equation.  Our  aim  In  this  chapter  is  mostly  to  Introduce 
principles;  the  methods  described  are  not  necessarily 
the  methods  to  be  used  in  everyday  practice.  More  advanced 
procedures  are  discussed  in  the  next  chapter. 

Section  3.2  The  Taylor  Series  Method  for  Initial 

Value  Problems 

As  a  preliminary  to  this  section,  we  first  describe 
the  reduction  of  higher-order  ordinary  differential 
equations  subject  to  initial  conditions,  to  sets  of  first 
order  equations. 

Although,  the  most  general  form  of  an  n!th.  order 
ordinary  differential  equation  is  that  of  (2.5*l)>  we  will 
consider  a  slightly  restricted  form,  namely 

Dny  (x)  =  f(Dn_1y(x),  Dn~2y(x) , . . . ,Dy(x) ,  y(x) ,  x) . 

(3.2.1.) 

In  other  words,  we  assume  that  we  can  solve 
(2.5.1)  formally  for  the  highest  derivative  present. 

Let  us  suppose  that  ( 3.2.1)  is  subject  to  the 
initial  conditions 

y(x  )  =  ri 
J  v  o  *o 

DyUo)  - 

Dny(x0)  =  nn 

at  the  point  x  =  xq . 


I 

' 
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We  write 

y0U)  =  y(x) 
yxU)  =  Dy(x) 
y2(x)  =  D2y(x) 

yn(x)  =  Dny(x)  =  f(yn_1>yn_2> • • • >y1>y0<x)  (3.2.3) 

Then,  differentiating  the  first  of  (3*2,3)#  we 

have 

Dy0  =  V  = y'  =  n- 

That  is 

yi=yl‘  <3-2-4) 

We  can  similarly  obtain 

y'i  =  y2 

y,2  =  y3 

•  * 

*  * 

*  » 

y'n-l  =  yn  =  f^yn-l’yn-2’ ' ‘ ' ,yo’x^ '  (3-2-5) 

The  initial  conditions  (3.2.2)  can  then  be  written 

yo(xo')  =  % 

yi(x0)  =  nx 

•  • 

=  1V  (3.2.6) 


' 


■ 
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It  is  then  clear  that  the  system 


y,1(x)  =  f1(x,y0,y1,y2. . . ,yn) 


(3.2.7) 


subject  to  initial  conditions  (3.2.6),  includes  (3.2.5)* 
and,  therefore,  we  are,  in  effect,  studying  the  solution 
of  (3.2.1)  and  (3.2.2). 

As  an  example  of  this  reduction,  consider  the 
non-linear  equation 


(3.2.8) 


X 


subject  to 


g(0)  =  0 

S'  (0)  =  T). 


(3.2.9) 


Writing 


yQU)  -  g(x) 
y1(x)=g1 (x) , 


(3.2.10) 


a 


■ 


■  - 


we  obtain  the  system 


y'0U) 
y' i(x) 


V*) 


1 


y2Q(x)/x2 


(3-2.11) 


subject  to 


y0(o)  =  0 
y1(o)  =  t) 


(3.2.12) 


Returning  our  attention  now  to  (3.2.7)*  we 
consider  a  typical  equation,  say  the  i 5  th. .  We  have 


y ' i (x)  =  f1(^yQ*y1* . . . *yn) .  (3.2.13) 


Differentiating  with  respect  to  x. 


<3f.  ,  <3f 
_ 1  + 

dx  c)y 


df.  .  ,  ,df.  , 

o  y  1+  + 


y",(x)  =  Xi  +  Xi  y',,  +  Xi  y'n  +  ...+Xi  .  (3.2.14) 


dyn  cFy 

o  1  n 


Assuming  sufficient  differentiability  of  f^  with 
respect  to  all  of  its  arguments,  we  can  go  on  to  obtain 
expressions  for  derivatives  of  y^(x)  of  all  orders.  To  do 
this,  we  must  proceed  recursively.  Using  (3.2.14),  we 
obtain  y".(x)  at  x  =  x  ,  for  all  values  of  i.  These  values 
we  substitute  into  the  expression  for  y?'(x),  obtaining 
y'"  (x  )  for  all  i,  and  so  forth,  for  as  many  terms  as  we 

i  O 
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need,  or  can  reasonably  obtain. 


We  will  write  y^(x)  as  a  Taylor  series-.  In 
particular,  for  x  =  xQ+h,  we  have 


y,  (x  +h)  =  y  (x  )  +  hy'  (x  )  +§Ty'.'(x^)  +.  .  . 


2! 


(3.2.15) 


Since  the  numbers  y^(xQ)  are  given  as  initial 
conditions,  and  the  values  of  all  required  derivatives  of 
y. (x  )  are  calculable  as  described  above,  we  can  obtain 
values  for  all  the  y. (x  +h) . 


Obviously,  we  do  not  intend  an  infinite  summation 
of  derivatives  in  (3-2.15).  Where  we  stop  depends  on  two 
things:  the  difficulty  of  evaluating  higher  derivatives; 

and,  the  necessity  of  obtaining  sufficient  convergence  with, 
the  chosen  value  of  h.. 


For  instance,  consider  the  function 

y  =  a  sin  x  +  bex,  (3-2.16) 

which  satisfies 

ylv  =  y.  (3.2.17) 

Suppose  the  initial  conditions  are  such  that  b  is  quite 
small.  Then  clearly,  the  solution  behaves  much  like  a 


■ 
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sine  function  for  small  x,  like  an  exponential  for  large 
x,  and  is  more  complex  for  intermediate  x.  For  a  given 
number  of  terms  in  the  Taylor  expansion,  it  is  clear  that 
a  different  h.  will  be  required  to  obtain  satisfactory 
convergence  in  each  region;  or,  conversely,  that  for  a  given 
fixed  h,  a  different  number  of  terms  in  the  Taylor  expansion 
will  be  required. 

In  practice,  we  often  discover  that  higher- 
order  derivatives  are  quite  difficult  to  evaluate,  and 
so  choose  h  small  enough  to  quarantee  convergence  to  the 
desired  number  fo  figures,  using  only  those  derivatives 
which  are  evaluable  with  reasonable  effort. 

Once  all  the  f . (x)  have  been  advanced  to  x=x  +h, 
we  can  regard  them  as  new  "initial  conditions,"  and 
advance  to  the  next  point,  usually  x  +2h.  We  continue  thus 
until  values  of  f^(x)  have  been  obtained  for  as  many  values 
of  x  as  desired. 

It  is  desirable  to  estimate  very  carefully  the 
amount  of  calculation  needed  to  advance  the  solution. 

Higher  derivatives  often  tend  to  have  many  terms,  and  thus 
are  costly  to  evaluate.  It  will  often  happen  that  a 
decrease  in  the  size  of  h.,  even  at  the  cost  of  more  points 
to  be  calculated,  will  give  accuracy  comparable  with 
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the  Inclusion  of  another  derivative  term,  and  at  less 
total  cost  in  effort  of  evaluation. 

Notice  that  there  is  no  difficulty  in  changing 
to  a  new  h.  in  mid-integration  if  the  convergence  either 
deteriorates  or  improves. 

The  Taylor  series  method,  as  described  above, 
has  advantages  and  disadvantages  over  other  methods  which 
will  be  discussed  in  this,  and  the  next,  chapter.  In  its 
favour,  it  is  applicable  to  a  very  wide  class  of  problems, 
including  non-linear  systems.  Also,  it  requires  no  special 
starting  procedure.  On  the  other  hand,  the  evaluation  of 
the  necessary  derivatives  may  often  be  difficult.  Further¬ 
more,  it  is  presently  unreasonable  to  expect  an  automatic 
digital  computer  program  to  evaluate  the  derivatives  of  an 
arbitrary  function,  given  only  its  formj  that  is,  without 
formulae  for  these  derivatives  being  worked  out  by  hand. 
However,  some  progress  has  been  made  in  this  area.  Also, 
for  some  special  classes  of  f^,  recursion  formulae  may  be 
available  for  the  derivatives. 

We  notice  also  that  the  Taylor  series  method  can 
be  applied  directly  to  second  and  higher  order  equations 

without  first  transforming  them  to  sets  of  first  order 

equations . 
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Section  3 o  Euler1 s  Method  for  Inital  Value  Problems 

Euler  [7l>  devised  a  method  for  the  approximate 
solution  of  systems  like  (3.2.7)  and  (3.2.6).  Before 
describing  it,  we  will  change  our  notation  to  make  it  more 
compact . 


We  will  write  y^  of  (3.2.6)  and  (3.2.7)  as 


components  of  a  vector 


X(x)  - 


^  y0(x)Ni 

yx(x) 


xyn<Xl  ^ 


(3.3.1) 


Hence,  we  can  write 


f1(x,yo,y1.  .  .  ,yn)  =  f\(x,2;) 


(3.3.2) 


Similarly,  we  write 


f0U>z) 
f\  (x,x) 


f  (x,x) 


fn(x,z) 

\  / 


(3.3.3) 


and  hence,  we  can  now  rewrite  equations  (3.2.7)  as  the 
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vector  equation 

y' (x)  =  f (x,y) . 


(3-3-4) 


If  we  also  write  the  initial  conditions  in 


vector  form,  as 


a  = 


'o 


1 


\  n  y 


(3.3.5) 


then  (3.2.6)  is  compactly  written 


(3.3.6) 


Operations  on  a  vector,  such  as  differentiation, 
differencing,  integration,  etc.,  are  to  he  interpreted  as 
acting  element  by  element  on  the  vector. 

Euler’s  method  can  now  be  obtained  by  at  least 
three  methods : 

(i)  Consider  the  Taylor  series  method  of  Section  3*2. 
Using  our  new  vector  notation,  equations  (3*2.15)  may  be 
written 

h2 

z(x0+b)  =  z(xQ)+b  y'(x)+  ^rx"(x)  +  ... 


(3.3-7) 


' 


U  -  ,  r 
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If  we  truncate  this  Taylor  series  after  the 
first  derivative  term,  and  substitute  (3*3 -4)  into  it, 


we  have 

x(xQ+h)^  y(xQ)  +hf(xo,y(xQ) ) . 

(3.3.8) 

(ii)  Alternatively,  we  may  integrate  (3.3-4) 

from 

x  to  x  +h 
o  o 


r  xo+h  xo+h 

J  Z  (t)dt  =  [  f(t,y(t))dt. 

xo  ^  o 

(3.3.9) 

That  is, 


xo+h. 

y(xo+h.)-y(xo)  =  fx  f(t,y(t))dt. 

(3.3.10) 

and,  approximating  the  integral  on  the  right  of  (3.3.10) 
by  the  product  of  the  interval  of  integration,  and  the 
value  of  the  integrand  at  the  left-hand  end  of  the  interval. 


z(xQ+h)  -  £(xo)  %  h.  f(xo,y(xQ)), 

(3.3.11) 

and  (3.3.8)  follows  immediately. 

(ill)  Finally,  we  use  the  approximation  (2.4.3)* 


h  D  *  A. 

^  rv  A 

That  is  D  ^  h  ^ 

(3.3.12) 

That  is 


(3.3.12) 


■ 
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and  thus. 


Equation  (3.3.8)  follows  immediately  from  the 
last  two  members  of  (3.3*13). 

The  formula  thus  derived  allows  us  to  calculate 
an  approximate  value  of  y;(xQ+h),  given  only  ;y(xQ)  and 
the  formula  for  f .  In  the  same  way  as  with  the  Taylor 
series  method,  we  can  advance  the  solution  a  step  at  a 
time  until  we  reach,  the  maximum  desired  value  of  x. 


We  exhibit  this  process  as  an  algorithm: 


1 


x«-x 


initialize  x 


o 


2 


h 


initialize  y_ 


3  PRINT;  x,y_ 

4  yf-y  +  h.  f(x,y) 


output  solution  point 


calculate  next  ordinate 


5 


X«-i x+h. 


calculate  next  abscissa 


6 


x  :x 


'max 


(<,>M3,7) 


test  for  maximum 
abscissa 


7 


STOP 


Oj'U' 
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The  Importance  of  this  algorithm  is  not  so 
much  practical,  since,  in  practice,  much  more  efficient 
algorithms  are  known,  but  is,  rather,  theoretical.  It 
is  the  basis  for  an  extremely  important  existence  theorem, 
which  we  will  discuss  in  the  next  section,  and  it  will 
also  serve  to  point  out  an  ubiquitous  difficulty  in  all 
initial  value  approximations  -  the  problem  of  stability  - 
to  be  taken  up  briefly  in  the  section  after  next. 

Section  3.4  Cauchy1 s  Existence  Theorem  for  Initial 

Value  Problems 

Subject  to  fairly  mild  restrictions  on  f(x,yj, 
Cauchy  [8]  proved  that  a  unique  solution  to  (3.3.4)  and 
(3.3.6)  exists,  and  that,  in  the  limit  as  h  — ►  0,  the  result 
of  Euler* s  method  converges  to  that  solution.  A  detailed 
account  of  this  proof  is  given  by  Henrici  [  9)  • 

As  a  simple  example,  we  consider  the  problem: 

y3  W  =  y(x)  (3.4.i) 

subject  to  the  initial  condition, 

y(o)  =  l.  (3. 4. 2) 

Let  us  suppose  that  we  use  Euler’s  method  for 
n  steps,  each,  of  length  h.,  arriving  finally  at  x  =  nh. 
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Then,  according  to  (3.3.8), 


y(xQ+h)  =  y(xQ)  +  hy(xQ)  =  (l+h)y(x0). 


(3.4.3) 


and,  in  general, 


y(o)  =  l, 
y(h)  =  (l+h), 
y(2h)  =  (l+h)2. 


y(nb)  =  (l+h) 


n 


(3.4.4) 


Now,  since  n  =  x/h,  the  last  of  (3.4.4)  gives 


y(x)  =  (i+h)x/h  =  (l+h)1/31 


x 


(3.4.5) 


Now,  as  h  -►  0,  by  the  definition  of  e, 

(l+h)1/11  —  e  =  2.71828...  , 


(3.4.6) 


and  so 


y(x)  —  ex. 


(3.4.7) 


which  is,  indeed,  the  exact  solution  to  the  given  problem. 
Thus,  the  convergence  of  Euler3 s  approximation  is  made 

plausible  at  least. 


Cauchy’s  theorem  can  be  stated  (and  proved)  in 
a  strong  form,  and  in  a  somewhat  weaker  form.  The  stronger 


, 


4 

,  (4+i )  i 

XI  r 


*  *  Ohi 


, 

s- 


08  0.:. 


,  's  -  •  (x  )y 


.  . '".r  c  .  ol'ii-  jjjs.Ec 


■■■  ■  &tfW9fiI08  '  ■ 


. 


-  44 


form,  as  proved  in  Henrici  [9],  follows. 

Consider  the  system 

z  M  =  >  (3.3.4) 

subject  to 

zU0)  =  a*  (3.3.6) 

Assume  that  f(x,y)  satisfies  the  following  two  conditions: 

(A)  f(x,y)  exists  and  is  continuous  for  x  e[a,b]  and 
y^,  the  components  of  y  in  (-00,00),  and,  (B)  there  exists 
some  constant  L,  called  the  Lipschitz  constant,  such  that 
for  x  e[a,b]  and  for  arbitarary  vectors  y  and  y*, 

||  f(x,y)  -  f(x,y*)  ||  <  L  ||  y-y*  ||  ,  (3.4.8) 

where,  for  any  vector,  say  v,  we  define  the  norm'1  of  v  as 

II V  II  =  |  vj-h  I  v2  |+  I  v3|+...  . 

(3.4.9) 

Under  assumptions  (A)  and  (B) ,  there  exists  a 
unique  vector  function  y(x)  which 

(i)  is  continuous  and  everywhere  differentiable  at 
least  once  for  x  €  [a,b], 

(ii)  which  satisfies  the  differential  system  (3.3.4) 


everywhere  for  x  e[a,b]. 
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and  (iii)  which  satisfies  the  initial  conditions  (3-3.6). 

This  theorem  is  stated  and  proved  in  a  weaker 
form,  namely  for  a  neighbourhood  only  of  x  =  a,  in 
Tricomi  [10.]  . 

Condition  (B)  is  called  a  Lipschitz  condition. 

It  is  a  condition  stronger  than  a  requirement  of  continuity 
alone,  but  not  as  strong  as  requiring  continuous  differen¬ 
tiability  . 


Leaving  the  vector  notation  behind  for  the 
moment,  and  considering  a  single  equation,  we  can  readily 
see  that  (B)  is  satisfied  by  any  continuously  differentiable 
f,  since,  by  the  mean  value  theorem. 


f(x,y)  -  f (x,y*) 

for  some  y  in  [y,y*] . 


6f(x,y) 

W 


(y-y*) 

(3-4.10) 


If  we  consider,  for  example  f(x,y)  =  jy|,  a 
function  which  is  not  continuously  differentiable,  we  find 
that  it  satisfies  condition  (B),  since  then 


ffcc.y)  -  f(x,y*)||=  I  (  ] y I  -  |y*|)|<|y-y* 


(3.4.11) 
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On  the  other  hand,  the  function  f(x,y)=|y|  '  ^ 
is  continuous,  but  does  not  satisfy  a  Lipschltz  condition. 
It  may  be  noted  that  for  functions  f  not  satisfying  a 
Lipschltz  condition,  a  solution  may  or  may  not  exist  or 
be  unique.  A  theorem  by  Peano  establishes  the  existence 
of  at  least  one  solution  (i.e.,  not  necessarily  unique) 
if  the  Lipschltz  condition  is  replaced  by  simple  continuity 

[17] . 


Section  3.5  The  Problem  of  Stability  of  Initial 

Value  Problems 

Cauchy’s  theorem  of  the  previous  section  might 
lead  us  to  believe  that  the  solution  of  initial  value 
problems  is  a  fairly  straightforward  procedure.  As  we 
shall  see  in  this  section,  such  a  belief  is  often  quite 
unfounded . 


Let  us  consider  the  second-order  differential 

equation 

y" ( x)  -  y(x)  =  0,  (3.5.1) 


subject  to 

y(o)  =  1,  and  y ' (0)  =  -1.  (3.5.2) 


The  exact  solution  of  this  system  can  readily 
be  found  to  be 

y(x)  =  e-x.  (3.5.3) 
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We  will  now  develop  a  numerical  method  of 
approximating  the  solution  to  this  system.  The  approach 
used  here  is  rather  atypical,  but  is  used  because  it  is 
somewhat  revealing. 

We  first  reduce  (3-5--1)  an^  (3-5.2)  as  described 
in  Section  3-2  to 

u  (x)  =  v(x) 

v  (x)  =  u(x) ,  (3-5-4) 

subject  to  u(o)  =  1  and  v(o)  =  -1, 
where  u(x)  =  y(x),  and  v(x)  =  y  (x) . 

The  Euler  approximation  gives  us 

u(x+h)  =  u(x)  +  bv(x)  (3-5-5) 

v(x+h)  =  v(x)  +  h.u(x) .  (3-5-6) 

Solving  (3-5-5)  for  v(x) ,  we  have 

v(x)  =  u(x+h)h-  u(x)  t  (3.5.7) 

A 

and ,  substituting  this  into  (3-5 •\)  > 

u(x+2h.)  -  u ( x+h. )  =  u(x+h.)  -u(x)  +  hu(x) 
h.  h.  \  J  • 

(3-5-8) 
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Hence,  we  obtain  the  second  order  difference 


equation 

u(x+2h)  -  2u(x+h.)  +  (l-h^)u(x)  =  0, 

(3.5.9) 


as  an  approximation  to  (3.5-1). 

Essentially,  (3.5*9)  allows  us  to  calculate 
the  next  step  in  the  solution,  given  the  previous  two. 

Now,  the  initial  conditions  (3.5.2)  give  us  the  value  of 
u(0).  Let  us  suppose  that,  before  we  start,  we  use  a 
Taylor  series  method,  or  some  other  method  to  obtain  an 
accurate  value  of  u(h),  or  simply,  suppose  that  u.(h)  is  given 
rather  than  u  (0) .  In  either  case,  let  us  use  as  the  value 
of  u(h.),  e  ,  the  exact  value  of  the  solution  to  the 
original  problem. 


The  results  of  a  numerical  calculation  using 
this  method  are  given,  for  three  values  of  h,  in  the  third 
column  of  Table  I,  along  with  corresponding  values  of  the 
exact  solution,  e-x,  in  the  second  column. 

The  values  obtained  are  drastically  in  error, 
expecially  for  h.  =  0.5.  Nor  h.  =  0.1,  the  approximation 
follows  the  correct  answer  fairly  well  for  a  while  but 
eventually  veers  off  and  grows  rapidly.  We  say  that  this 
approximation  is  "unstable." 
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TABLE  I 

SOLUTION  OF 

u(x+2h.)  =  2u(x+h.)  -  (l-h2)u(x) 


X 

e~x 

u(x)  using 
u(b)  = 

u(x)  using 
u(b)  =  (l-h) 

b 

=  0.5 

0.0 

1.Q0Q0Q 

l.QOpOQ 

1.00000 

0.5 

. 60653 

.60653 

.50000 

1.0 

.36788 

.46306 

.25000 

1.5 

.22313 

.47123 

.12500 

2.0 

.13534 

.59516 

.06250 

2.5 

.08209 

.83689 

.03125 

3.0 

.04979 

1.22741 

.01563 

3.5 

.03020 

1.82716 

.00781 

4.0 

.01832 

2.73376 

.00391 

b  -  0.25 

0.0 

1 . 00000 

1 . 00000 

1 . 00000 

0.5 

.60653 

.62010 

.56250 

1.0 

.36788 

.43881 

.31641 

1.5 

.22313 

. 38746 

.17798 

2.0 

.13534 

.43768 

.10011 

2.5 

.08209 

.58953 

.05631 

3.0 

.04979 

.86807 

,03168 

3.5 

.03020 

1.32651 

.01782 

4.0 

.01832 

2.05588 

.01002 

i — 1 

o 

ii 

0.0 

1.00000 

1.00000 

1.00000 

0.5 

.60653 

.61516 

.59049 

1.0 

.36788 

. 40298 

.34868 

1.5 

.22313 

.30194 

.20589 

2.0 

.13534 

.28135 

.12158 

2.5 

.08209 

.33210 

.07179 

3.0 

•04979 

.46339 

.04239 

3.5 

.03020 

.70411 

.02503 

4.0 

.01832 

1.10906 

.01478 
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To  see  why  this  is  the  case  we  will  solve  the 
difference  equation  (3.5*9)  analytically.  To  do  so,  let 
us  set 

x  =  nh.  and  u(x)  =  u(nh)  =  un .  (3-5.10) 

Then  we  may  write  (3.5*9) 

u  ,  n  -  2u  ,  ,  +  (l-h2)u  =0.  (3.5.11) 

n+2  n+1  v  '  n  w  ^  y 


As  a  trial  solution,  let  us  use 


pn 

u  =  C  , 
n 


for  some  constant,  C. 


Then 

Cn+2 -2Cn+1+ ( 1 -h2 ) Cn=0 

and  so 

C2-2C+(l-h2)  -  0, 
of  which  the  roots  are 

C  -  1+h. 


Hence,  we  shall  write 

u  =  A(l+h)n+B(l-h)n 


(3.5.12) 


(3.5.13) 


(3.5.14) 


(3.5.15) 


(3.5.16) 


To  determine  the  coefficients  A  and  B,  we  will 
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use  the  initial  values  u(o)  =  1,  u(h)  =  e  h,  to  obtain; 
for  n  =  0, 

1  =  A+B  (3.5.17) 


and  for  n  =  1, 

e-h  =  A(l+h)+B(l-h) . 


(3.5.18) 


From  the  first  of  these, 
A  =  1 -B, 


and  therefore, 

e"h  =  (l-B)(l+h)+B(l-h) . 


(3.5.19) 


(3.5.20) 


Hence , 


B 


h.+  (l-e~h  ) 
2h  ' 


(3.5.21) 


and  therefore,  , 

A  .  h-(l-e-h)  _ 

2h 


(3.5.22) 


The  source  of  the  difficulty  is  now  clear.  The 
differential  equation  (3.5-1)  has  both  an  increasing  and 
a  decreasing  solution,  but  the  Initial  conditions  specify 
a  coefficient  of  zero  for  the  increasing  solution.  The 
difference  equation  (3-5.12)  also  has  increasing  and 
decreasing  solutions,  (l+h)n  and  (l-h.)n,  but  the  initial 
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conditions  that  we  used  allow  a  component  of  the  increasing 
solution  to  remain.  This  component  quite  quickly  dominates 
the  decreasing  solution,  with  the  results  which  we  have 
observed . 

We  can  see  that,  in  a  more  general  situation, 
the  difference  equation  might  have  several  increasing,  or 
several  decreasing  solutions,  and  that  a  component  of  an 
unwanted  solution  might  creep  in  and  swamp  the  desired 
solution . 

It  can  be  seen  that  if  we  had  used  u(h)  =  1-h. 

— h 

instead  of  e  ,  that  we  would  have  obtained  A  =  0  and 
B  =  1.  In  other  words,  the  increasing  solution  would  have 
been  excluded.  The  fourth  column  of  Table  I  lists  solutions 
of  (3-5.9)  using  u(h)  =  1-h.  We  observe  a  very  much 
improved  behaviour,  although  the  approximation  is  still 
far  from  excellent. 

It  is  interesting  to  note  that  had  we,  in  obtaining 

our  value  of  u(h.)  been  less  fastidious  and  used  (3.5.5) 

.t"-) 

instead  of  the  true  value  u(h)  =  e  ,  we  would  have  obtained 
u(h)  =  1-h. 


We  have  now  illustrated  the  fact  that  the  exact 
initial  conditions  of  the  differential  equation  may  often 
(in  fact,  almost  always)  not  be  the  most  appropriate 
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conditions  to  use  in  a  difference  approximation  to  that 
differential  equation. 

Also  worth  noting  is  that,  even  using  u(h)  =  1-h, 
we  expect  round-off  error  in  our  arithmetic  to  introduce 
a  small  component  of  the  increasing  solution,  and  that, 
eventually,  the  decreasing  solution  would  become  dominated. 

However,  this  situation  is  far  less  serious  than  using 

— h 

u(h)  =  e-  ,  as  Table  I  shows. 

The  topic  of  stability  is  a  large  one,  as  we  shall 
return  to  it  briefly  in  later  sections. 

Section  3.6  Simple  Approximations  for  Boundary 

Value  Problems 

In  this  section,  we  will  consider  linear  dif¬ 
ferential  equations  subject  to  boundary  conditions. 

Existence  theorems  are  available  for  these  systems  [9l>[l0]. 
For  instance  [9],  for  the  problem 

y"  =  f (x,y) ,  (3.6.1) 

subject  to 

y (a)  =  a,  y(b)  =  p,  (3.6.2) 

there  exists  a  unique  solution  y,  provided  that  df (x,y ) 

dy 

is  continuous  and  non-negative  for  x  e[a,b],  and  y  e(-00,00). 
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Although  most  of  the  results  of  this  section 
can  readily  be  extended  to  higher-order  equations,  we 
will,  for  simplicity,  consider  only  the  following  second- 
order  equation: 

y"(x)  +  f(x)y'(x)  +  g(x)y(x)  =  k(x). 

(3.6.3) 

We  will  first  suppose  that  we  have  boundary 
conditions  of  the  form  (3.6.2).  We  will  attempt  to 
calculate  approximate  values  y^  of  y(x)  at  the  n+1  points 

x  =  a,  a+h,  a+2h, . . . ,  b-h,  b;  (3.6.4) 

we  write 

yi  *  y(a+ih) .  (3.6.5) 


Let  us  begin  by  replacing  the  derivatives  in 
(3.6.3)  by  the  difference  expansions  (2.3.11)  and  (2.3.15); 
each  truncated  to  a  single  term;  thus 

h2D2  %  52 ,  (3.6.7) 

and 


h.D  ^  [_l6  . 


(3.6.8) 
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Then,  we  will  have 


h2D2yi  %  b2Yi  =  y1+i-2y1+y1-1. 


(3.6.10) 


and 


hDyi  *  n6y±  =  |(y1+1-y1.1); 


(3.6.11) 


which,  when  substituted  into  (3.6.3) >  yield 


(yi+i  2yi+yi-i^  +  -2i(yi+i-yi_i)  +  h2siyi  =  h2k±. 


i.e . , 

hf.  p  hf 

yi+i(l+— )  +  y^h  g.-2)  +  y^d-  —  )=  h  k. . 

(3.6.12) 

Writing  down  (3.6.12)  for  successive  values  of 
i,  and  keeping  the  known  boundary  values  in  mind,  we  have 


-  ( 2-h2g1  )yx+(  1+^-h.f  1  )y2 


:hdk1-a(l-^hf1) 


( 1  -|hf 2 ) y  1  -  ( 2  -h2^  )  y2+  ( l+|hf  2 )  y^ 


p 

=h.^k. 


( 1 -|b.f 3 ) y2 - ( 2 -h2g3 ) y3+ ( l+|hf 3 ) y^  =h2k 


3 


(3.6.13) 


(Cl  £.-■) 
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( 1  -4hf  n)y  0-(2-h‘"g  „ ) y  ^(l+^-hf  0)y  ,  =h^k  „ 

N  2  n-2/‘yn-3  &n-2 '  Jn-2  v  2  n-2/,yn-l  n-2 

(l-feif  )y  0-(2-h2g  ,  )y  ,  =b2k  ,-p(l+fai.f  n), 

v  2  n-l/l/n-2  v  &n-l/Jn-l  n-1  2  n-l'* 

a  set  of  (n-l)  linear  equations  in  the  (n-l)  unknowns 


yi'y2' • . • ^yn_x- 


Let  us  write  (3-6.13)  in  matrix-vector  form 


A  =  d , 


(3-6.14) 


where  A  is  the  (n-l)  hy  (n-l)  square  matrix 


-(2-hig1) 


hf 

(1— 2s) 
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-(2-h2g2) 
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0 


hfo 

d+  -5-) 


hf. 


-(2-h2g3)  (1+  — 2^-) 


h.f  q  o  hf  q 
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y  is  the  column  vector 


/ 


y- 

y. 


N 


kyn-y 


and  d  is  the  column  vector 


/ 


P  hf 

h  k1-a(l-  ) 


h2k. 


2 

h^k. 


h^k 


\ 


\ 


n-2 


o  hf  n 

h2kn.i-P(i+  -§=i> 


Since  the  matrix  A  is  triple-diagonal,  (3.6.14) 
may  most  readily  be  solved  by  the  LU  factorization  method 
[ 11 ] .  When  this  method  is  programmed  for  a  digital 
computer,  we  note  that  storage  requirements  go  up  only 
linearly  (rather  than  quadratically )  and  the  computing 
time  also  goes  up  linearly  (rather  than  cubically)  with 
the  number  of  equations  to  be  solved. 
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Suppose  that,  instead  of  boundary  conditions  of 
the  form  (3.6.2),  we  are  given 


y  (a)  +  7y( a)  =  a 
in  place  of  y(a)  =  a. 


(3.6.15) 


Then,  using  (3-6.12)  with  i  =  0,  we  have 

hf  0  hf  n 

(l - 2^)y_1-(2-h2go)yQ+(l+  )y1  =  h2kQ 

(3.6.16) 

and  the  finite  difference  replacement  for  (3-6.15)  is 


y1-y_1+2h.yyo  =  2ha. 


(3.6.17) 


we  obtain 
hf 


Eliminating  y_-^  between  (3.6.16)  and  (3. 6.17), 


hf 


[(l-  —rp~) 2hy- ( 2-h2g^ ) ]  y^+2yn  =  h2k^+2ha(l - 5—). 

(3.6.18) 


This  becomes  an  additional  first  equation  of  the 
set  (3.6.13),  and  the  second  equation  (replacing  the  old 

first  equation  )  is  clearly 

hf  hf 

(1-  -gi)yo-(2-h  g1)y1+(l+  )y2  =  h  . 


(3.6.19) 
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The  remaining  members  of  (3. 6.13)  are  unaltered 


unless  the  boundary  condition  at  x  =  b  also  contains  a 
derivative,  in  which  case  (3.6.13)  acquires  an  additional 
last  equation,  and  the  old  last  equation  is  modified  in 
the  obvious  way. 

Finally,  we  will  consider  briefly  a  class  of  eigenvalu 
differential  equations .  If  we  are  given 


y"(x)+f(x)y  (x)+g(x)y(x)  =  Xy(x),  (3-6.20) 


subject  to 


y(a)  =  y ( to )  =  0 


(3.6.21) 


then  we  can  perform  the  analysis  of  the  beginning  of  this 


section.  When  we  do  so,  we  obtain  a  system  of  equations 

like  (3.6.13)  or  (3.6.14)  except  for  the  right  hand  side 

2 

d_.  Clearly,  d_  is  replaced  by  h  Xy. 


Thus,  we  have  to  solve  the  problem 

2. 

A  =  h  Xy_  , 


(3-6.22) 


or,  setting 


(3.6.23) 


we  have  to  solve 


A  £  =  U-2. 


(3-6.24) 
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Equation  (3.6.24)  is  the  standard  algebraic 
eigenproblem,  and  may  be  solved  by  standard  techniques 
which  need  not  concern  us  here. 


Section  3.7  Classification  of  Second-Order  Partial 

Differential  Equations 

We  will  consider  the  general,  second-order 
quasi-linear  partial  differential  equation: 


-n2  >2 

o  u  ,  ,  o  u 

- Q  +  t) 

6x2 


dxdy 


+  c 


02u 

6y2 


-  e. 


(3.7.1) 


6u  6u 


where  a,b,c,  and  e  are  real  functions  of  u,^-, x,  and  y, 
but  not  of  the  second  derivatives  (hence  the  classification 
"quasi-linear" ) . 


The  following  discussion  could  equally  well  be 
carried  out  with  respect  to  a  pair  of  first-order  equations 
equivalent  to  (3«7.l)>  instead  of  the  single  second-order 
equation . 


A  standard  notation  to  compact  (3.7-1)  is 


6u  6u 

p  =  5P  q  =  3y'  r  = 


whence  (3*7 • l)  becomes 


a!u 

Sx2 


^2 
o  u 

’  s  =  = 


32u 


6y 


2  ' 


(3-7.2) 


=  e . 


ar  +  bs  +  ct 


(3.7.3) 
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The  classification  of  (3-7 »3)  is  motivated  by 
a  simple  attempt  to  solve  it.  Let  us  suppose  that  we  are 
given  starting  data  in  the  form  of  values  of  u,p,  and  q 
along  some  curve,  say  L,  in  the  (x,y)  plane  described  by 
a  parameter  i.  That  is  to  say,  we  are  given  u,  p,  and  q 
as 

u  =  u(i),  P  =  p(iO>  and  q.  =  q(j&)> 

(3.7.4) 


on  the  curve  L  defined  by 


x  =  xU)  and  y  =  y(i).  (3-7.5) 


Our  attempted  solution  is  by  means  of  a 
generalization  of  the  Taylor  series  method  of  Section 
3-2.  Our  starting  conditions  give  us  the  value  of  u  and 
its  two  first  derivatives  at  any  point  on  L.  Using  this 
information,  and  the  differential  equation  (3*7.3)*  let 
us  attempt  to  find  values  for  the  three  second  derivatives, 
four  third  derivatives  and  so  on,  of  u. 

By  simple  calculus,  we  can  write 

dp  =  ||dx  +  ^  dy,  (3-7-6) 

or 


dp  =  rdx  +  sdy. 


(3.7.7) 
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If  we  differentiate  along  L,  i.e.,  with  respect 


to  & ,  then 


_  r  dx 

di  "  di  °  an  * 


(3-7.6) 


Similarly,  by  differentiating  q  along  L, 


s  6x 
di  t 


d n 

an 


(3.7.9) 


Thus,  (3.7.3)..  (3.7.6),  and  (3- 7. 9)  are  three 
equations  in  the  three  unknowns  r,  s,  and  t  which  we 
seek.  We  write  these  three  equations  in  matrix  form  as 


/  a 

b 

c  > 

frNi 

e 

dx 

an 

dy 

an 

0 

s 

= 

dp 

an 

(3.7.10) 

0 

\ 

dx 

an 

dy 

an  ) 

vV 

dq 

y 

A 

f  = 

s  • 

(3.7.11) 

If  we  suppose  that  we  have  solved  (3-7.ll)  for 
r,  s,  and  t,  then  we  can  differentiate  (3 -7.il)  with 
respect,  say,  to  x,  to  obtain 


a  +  14  f 


ox 


ox 


dx 


(3.7.12) 


and  we  then  solve  (3.8.12)  for  8f  . 


If  we  have  already  obtained  A  1  in  solving 
(3-7.ll) ,  then  we  may  use  it  again,  and  write 


8f  _  A-1  dg. 

_ 


8a  f 


-  8x  8x 


(3-7.14) 


The  elements  of  _=  are,  obviously, 

8x 

Sr  _  A,  Ss  _  o3u  _  St  =  —  u  p  ,  (3*7.15) 

37  “  Sx3  37  "  Sx2Sy  ;  37  3x3y 

three  of  the  required  four  third  derivatives  of  u. 

The  fourth,  can  be  found  by  differentiating 
(3.7.11)  with  respect  to  y,  and  proceeding  as  above. 

In  this  fashion  we  can,  at  least  in  principle, 
go  on  to  obtain  derivatives  of  all  orders,  and  carry  out 
our  Taylor  series  scheme. 


Unfortunately,  there  is  an  implicit  assumption 
involved  here,  namely  that  the  matrix  A  is  non-singular. 
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Since  the  manifold  L  has  so  far  been  arbitrary,  and  since 
a,  b,  and  c  are  arbitrary,  there  is  no  reason  why  A 
should  not  be  singular. 


If,  indeed,  A  is  singular,  then 


a|  =  a(^y  -  -  0 

-i  '■cU'  0df  df  'an1  ~  ’ 


(3-7.16) 


which  we  rewrite  as 

2 


-»S+-o, 


(3.7.17) 


or,  finally,  as 


dy 

dx 


^b^-4ac 

2a 


(3.7.18) 


Clearly,  if  A  =  b  ~4ae  is  positive,  then  (3.7*18) 
defines,  at  any  given  point  on  L,  two  directions  along 
which  the  Taylor  series  scheme  is  certain  to  fail.  If 
A  =  0,  there  is  only  one  such  direction,  and  if  A<0, 
there  are  no  such,  directions  in  the  real  domain. 


According  as  A>0,  A=0,  or  A<0,  we  term  the 
differential  equation  (3*7.3)  "hyperbolic,"  "parabolic," 
or  "elliptic,"  respectively  [12]. 

It  is  interesting  to  notice  that,  since,  in 
general,  a,  b,  and  c  may  depend  upon  u,  the  character 
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(i.e.,  hyperbolic,  etc.)  of  a  given  equation  may  change 
from  region  to  region. 

Either  of  (3. 7*17)  or  (3*7.18)  is  called  the 
"characteristic  equation"  of  (3*7*3)  and  the  solutions 
are  called  "characteristics"  or  "jump  lines."  The  reason 
for  the  latter  term  is  clear,  since  the  jump  lines  are 
lines  along  which  the  Taylor  series  for  u  may  fail  to 
exist,  and  hence,  along  which  u  may  be  discontinuous. 

This  property  is  elucidated  in  some  detail  in  [13.1*  page 
4l6 . 


It  can  be  shown  that  the  type  of  a  partial 
differential  equation  determines  to  a  very  great  extent 
the  type  of  boundary  conditions  which,  are  required  to 
make  the  equation  solvable  in  a  realistic  manner.  A 
discussion  of  this  topic  may  be  found  in  [12],  or,  at 
greater  length,  in  [13] • 

We  will  discuss  typical  methods  for  solving 
systems  of  the  three  types,  in  the  succeeding  three  sections. 

Section  3.8  The  Method  of  Characteristics  for 

Hyperbolic  Systems 

The  characteristic  curves  described  in  the 
previous  section  prevented  our  easy  exploitation  of  the 
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Taylor  series  approach  to  a  solution.  In  this  section, 
we  make  use  of  these  same  characteristics  to  obtain  an 
approximate  solution  to  a  hyperbolic  equation  of  the  form 
(3.7.3).  The  notation  and  exposition  of  this  section 
follow  that  of  [ll] . 


Referring  to  Figure  I,  let  us  suppose  that  we 


are  given  the  values  of  u,  p,  and  q  on  a  segment,  AB,  of 
some  non-characteristic  curve.  Considering  any  two  points 
P  and  Q  on  AB  (which  are  fairly  close  together),  then 
two  characteristics  pass  through  each  of  P  and  Q.  These 
curves  satisfy 


b^-4ac 


(3.7.18) 


2a 

and  we  shall  distinguish  the  two  cases  by  writing  them  as 


(3.8.1) 


where 


f  =  b*K/b2-4ac 


(3-8.2) 


Returning  now  to  (3*7*ll)>  since  we  are  on  a 


characteristic,  A 


=  0.  Since  we  are  seeking  an  analytic 
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solution,  it  is  clear  that  for  such  a  solution  to  be  possible, 


the 

quantities  r,  s 

,  and 

t  must  none- 

-the-less  exist.  That 

is 

if  we  were  to  attempt 

to  solve 

(3.7.11) 

by  Cramer’s  rule. 

then  we  would  have 

for  r 

9 

e 

b 

c 

d£ 

d^ 

0 

di 

d  i 

dq, 

dx 

dZ 

d  l 

d  i 

d  l 

r  = 

(3.8.3 

|  A  | 

Hence,  if 

r  is 

to  exist 

,  we  must 

have 

e 

b 

c 

dp 
d  i 

dZ 
d  i 

0 

= 

0,  also. 

dq 

dx 

dZ 

d  i 

dZ 

(3.8.4) 

Similarly 

,  in 

order  that  s 

and  t 

exist , 

a 

e 

c  1 

dx 

0 

<u 

d  i 

= 

0, 

0 

dq 

dx 

(H 

dZ 

(3.8.5) 
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and 


a 

b 

e 

dx 

dZ 

dp 

d  l 

0 

dx 

dq 

d^ 

di 

=  0.  (3-8.6) 


The  previous  three  equations  are  equivalent 
on  a  characteristic  and  we  arbitrarily  choose  (3-8.5) , 
to  give 


dx  dy  „  dp  dy  dx  dq  ~ 

61  61  a  61  61  0  61  61  ~0, 


(3-8.7) 


which  we  rewrite  as 

e  dy  -  a  dp  ^  -  c  d  q  =  0.  (3-8.8) 

Let  the  characteristic  of  the  f  type  through 
P  meet  the  characteristic  of  the  g  type  through  Q  at 
some  point  S.  Then,  to  a  first  approximation,  we  regard 
PS  as  a  straight  line  with  slope  f p .  Hence,  approximately, 

ys  -  yp  =  fp(xs_xp)>  (3-8.9) 

and,  by  similarly  regarding  QS  as  straight,  with  slope 

i’ 

yS  -  yQ  =  •  (3.8.10) 
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Equations  (3.8.9)  and  (3.8.10)  are  two  linear 
equations  which  we  may  readily  solve  for  x^  and  y^,  since 
all  the  other  quantities  are  known. 


Hence,  using  (3.8.8)  along  the  lines  PS  and  QS, 
and  using  our  initial  estimates  of  and  y^,  we  can 
write 


ep(ys  Yp)  (p g ~Pp ) Pp -Cp ( q^ ~qp )  —  0, 

(3.8.11) 


and 


eQ(ys-yQ)-aQ(ps-pQ)sQ-cQ(qs-qQ)=  0. 


(3.8.12) 


Again,  (3.8.11)  and  (3.8,12)  are  two  linear 
equations  which  we  can  solve  for  and  q^ .  Finally,  we 


may  obtain  u^  approximately  from 
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where  by 
Hence , 
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ou 

and 

we  intend  mean  values  along  PS. 
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2 
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(3.8.14) 


Having  now  obtained  an  approximate  value  of 
u  at  S,  we  can  improve  it  by  the  following  process: 
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(l)  Substitute  the  approximate  values  of  Ug,  Pg, 
and  qg  Into  (3.8.2)  to  obtain  approximate  values  of 
fg  and  gg . 

(ii)  In  place  of  (3-8.9)  and  (3-8.10),  use 

ys-yp  =  (  Sg— ■ - ) (xs-xp) ,  (3.8.16) 

and 

yS‘yQ  =  (^S-)(xS-xQ)’  (3-8.17) 

to  obtain  improved  values  of  Xg  and  yg. 

(iii)  In  place  of  (3-8.ll)  and  (3-8.12),  use 


0  Q1 


(ys-yp) 


3-p  3-^ 


( Pg  pp ) 


rfp+fg' 


rCp+Cg 


( q<^  Q.P )  ~  0 


(3-8.18) 


and 


raQ+aS 


2  I^S-V 


r  fQ+fS 


CQ+C'S 


(qs-qQ)=  0, 


(3.8.19) 


to  obtain  improved  values  of  Pg  and  qg . 
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(iv)  Use  (3.3.14)  to  obtain  an  Improved  u0  • 

Steps  (l)  to  (iv)  are  to  be  repeated  until  two 
successive  values  of  ug  agree  to  the  desired  accuracy. 

In  a  similar  fashion,  using  the  points  Q,  and  R 
in  Figure  I,  we  can  obtain  u^  and  its  derivatives.  Then 
building  out  from  S  and  T,  we  obtain  the  solution  at  V. 

Finally,  we  can  build  up  a  complete  solution  in 
the  region  ABC.  It  is  clear  from  our  method  of  solution 
that  we  can  go  no  further,  because  points  outside  of  ABC 
depend  on  values  of  u,  p,  and  q  on  the  continuation  of  AB. 
Furthermore,  it  is  clear  that  the  solution  Inside  of  ABC 
is  entirely  independent  of  such,  values. 

Section  3.9  Solution  of  Parabolic  Systems 

In  this  section,  we  will  only  illustrate  the 
solution  of  parabolic  equations  by  consideration  of  a 
specific  parabolic  equation,  the  dimensionless  form  of 
the  heat  equation: 

_  df  .  (3.9.1) 

3x2  “  *£ 

As  subsidiary  conditions,  f  is  specified  on 
three  sides  of  a  rectangle,  say  t  =  0,  and  x  +  +1 . 


he..  OTpfiil  •-$  n  ivtrlc  o.:!  (. ^  .v .  8 . )  -laU 

...  ;iyd*8 


. 


-  7.3 


Probably  the  simplest  procedure  for  parabolic 
equations  Is  to  construct  a  finite  difference  approximation 
to  the  second  derivative  term,  thereby  generating  a  set 
of  initial  value  equations. 


Dividing  the  range  [-1,1]  Into  equal  intervals 
of  width  h,  with  tabular  points  x  ,x^...,xn,  and  denoting 
f  evaluated  as  a  function  of  t  for  x  =  x^  by  f^(t),  we 
can  replace  (3.9-1)  by  the  equations 

h2  dfi(t)  =  f.  i , ;  1=1, 2,..., n-1 
dt  liii+i 

(3.9.2) 

2 

where  the  right  hand  side  is  the  simple  5  approximation 
to  the  second  derivative. 


Writing  down  (3-9.2)  for  all  values  of  i  and 
including  the  boundary  conditions,  we  have 


df 

2  1  +  2f  -f 

h2  _1  2fx  f2 


2  d^2 

h  — —  -  f  +2f  -f 
dt  1  d  2  3 


f 

o 


=  0 


h 


2  df3 
dt 


-  f2+2f3-f4 


=  0  (3-9-3) 


h. 


2  dfn-2 
dt 


■f  0+2f  0-f  i 

n-3  n-2  n-1 


=  0 


h. 


o  df  n 
2  n-1 

dt 


-f  0+2f  , 

n-2  n-1 


=  f 


'n 
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The  subsidiary  conditions  provide  all  the 


f  at  t=0,  and  f  and  f  for  all  t.  Hence  (3. 9 -3)  is 
r  on 

a  system  of  the  type  of  (3*9.7)  and  (3*2.6)  and  is 
solvable  by  the  techniques  of  Sections  3*2  and  3*3>  or 
by  the  more  sophlstocated  techniques  described  in  the 
next  chapter.  It  is  also  subject,  of  course,  to  the 
same  stability  problems  that  we  saw  in  Section  3*5* 

Section  3*10  Solution  of  Elliptic  Systems 

Boundary  conditions  for  elliptic  equations  must 
be  given  on  a  closed  contour.  The  usual  method  of 
solution  is  to  divide  the  region  inside  the  contour  into 
a  grid  of  square  or  rectangular  cells,  replace  the 
differential  equation  by  finite  difference  equations,  and 
solve  the  resulting  algebraic  equations. 

By  way  of  illustration,  we  will  consider  the 
Laplace  equation. 


(3.10.1) 


with,  values  of  f  given  on  a  rectangular  boundary. 

Referring  to  a  particular  point  (xo,yQ)  of 
the  grid,  we  may  write 


o 


(3.10.2) 


, 


M  •  :  :  o  )0$  G  Jt  , 

'  •-#  M'  f  i:  '•  3V  -  j-  ■  9l;  ./J  [•; 


where  the  subscript  x  on  5  implies  differencing  in 
the  x-direction.  A  similar  expression  can  be  written 
for  the  derivative  with  respect  to  y,  whence  we  will 
write 


2 

5  i 
x  o 


=  (fn  -2f  +f  ,  ) 

'  1  o  -l'x 


o 

5  “f 

y  o 


=  (f,  -2f  +f  )  . 
'  1  o  -1  y 


(3.10.3) 


We  thus  replace  (3.10.1)  hy 


(ff2Vf-l>x  +  <§>  (fr2fo+f-l  >y  -  °-  (3.10.4) 


which,  in  the  case  5x=6y,  becomes 


(f  +f  i )  +(f,+f  a )  -4f  =0 
VX1  -lyx  v  1  -1  y  o 


(3.10.5) 


Thus,  we  can  set  up  an  equation  like  (3-10.5) 
at  every  point  of  the  grid.  If  we  have  a  hundred 
subdivisions  in  the  x  and  the  y  direction,  then  we  will 
have  ten  thousand  simultaneous  linear  equations.  In  as 
many  unknowns,  to  solve.  Practical  techniques  for  solving 
such  huge  sets  of  equations  on  a  digital  computer  are 
described  by  Varga  [l4]  . 


A  technique  for  handling  non-rec tangular 


' 
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boundary  contours  is  described  by  Murray  [15]. 
Detailed  discussions  of  these  and  many  other  problems 
associated  with  partial  differential  equations  may  be 
found  in  [12],  [13]  >  and  [19]  • 
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CHAPTER  IV 

ADVANCED  NUMERICAL  METHODS  FOR  ORDINARY 
DIFFERENTIAL  EQUATIONS 


In  this  chapter  we  develop  and  discuss  practical 
methods  for  the  approximate  solutions  of  ordinary 
differential  equations. 

General  multi-step  methods  are  discussed, 
including  forward  integration  and  predictor-corrector 
formulas.  The  method  of  Milne  is  developed  in  detail,  and 
an  algorithm  for  it  is  presented. 

The  stability  of  multi-step  methods  is  analyzed 
in  some  detail.  Milne’s  method  is  shown  to  be  unstable,  but 
Hamming’s  modification  to  Milne’s  method  is  presented  and 
shown  to  be  conditionally  stable. 

Runge-Kutta  methods  are  introduced,  and  a  second 
order  formula  is  developed.  An  algorithm  for  the  Runge- 
Kutta-Gill  fourth  order  process  is  given. 

Difference  correction  methods  and  a  method  of 
linear  combinations  of  initial  value  solutions  for  boundary 
value  problems  are  outlined. 


Section  4.1  Introduction  to  Multi-step  Methods 

We  consider  in  this  section  a  large  and  very 
general  class  of  numerical  formulas  for  initial  value 
problems  -  the  multi-step  methods.  The  problem  to  be  solved 
is  again 

Z  M  =  ,  (3-3.4) 


subject  to 


s.(x0)  =  a- 


(3.3.6) 


'■  -  - ^  :  r:;  0iiAt .J  a 

■  -  '  _  :■  L .  .. 


t  '  J  '  ’  ov.  -"6  ,.,.0,:  ■  ,  ■; 

.  v::  a  t;.;  I  .  1  xTl  _• 

'  '  .  . 

-  '  '  '  '  •  '  •  '  •  '  -  ■  '  ■  -i  j'-'-  ■  '■  •■■.;  :■•  r... . ,  • 

.  ■ 

■  •'■■■  ■  >  '  :  :  X  :  i  ■  :  ~  ,,:b 

:  ■  lir  .  .  :h  ; ,  ■„  r,  w:0?, 

c 


..  ,  '■  .  :>  -  <  7  J  f  to*?. 

■  .  ■  •  - .  . .  ,  ■ .  - 

»  ■'  ?  ■■■'■.:  .  :  :  j:..  V"  :■  Bft 

T'J 7  X<-  ■;?  ;H  r  ,  J;.  ,;V;  ■; 

. 


*  ■  _ ) 


-  78  - 

For  simplicity  of  presentation,  we  will  consider 
only  a  single  first  order  equation,  Instead  of  a  set  of 
such  equations.  All  the  results  of  this  section  can 
be  extended  in  a  trivial  fashion  to  sets,  however.  Thus, 
we  consider 

y* (x)  -  f (x,y) ,  (4.1.1) 

subject  to 

y(xo)  =-yQ.  (4.1.2) 


The  name  multi-step  arises  from  the  form  of 

the  formulas  -  the  approximate  value  y^  of  y(xn)=y(xo+nh.) 

! 

is  expressed  as  a  linear  combination  of  y  and  y^  for 
several  values  of  i<n.  More  specifically,  we  can  consider 
two  types  of  multi-step  formula,  the  "forward  integration" 
and  the  "iterative"  types: 

(i)  In  forward  integration  formulas,  yn+1  is  expressed 

i 

as  a  linear  combination  of  y^  and  y^  for  several  values 
of  i<n; 

! 

(ii)  In  iterative  formulas,  yR+1  is  included  in  the 
linear  combination. 


In  general  then,  we  will  write 
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y 


n+1  Jrrvn4 an-l^n-l *  *  *  '  tan-p^n-p 


+h(b  , ,  y  , t +b  y  + . . . +b  y  ) 
n+l^n+1  n  n  n-p^n-p' 


+E  , 
n* 


(4.1.3) 


where  bn+-^=0  for  a  forward  integration  formula. 
represents  the  error  at  the  (n+l)th.  step,  supposing  that 
all  the  y^  and  y^  on  the  right  in  (4.1.3)  are  known  exactly 
(4.1.3)  is  to  be  regarded  as  an  equation  for  determining 
yn+1  explicitly  or  implicitly. 

The  error,  E'n,  arises  from  two  sources: 

(i)  Sincethe  right-hand  side  of  (4.1.3)  is  not  an 
infinite  sum  of  terms,  it  cannot  in  general  represent  yn+-]_ 
exactly.  That  is,  (4.1.3)  may  be  thought  of  as  being 
obtained  by  truncating  some  infinite  series.  This  type  of 
error  is  therefore  called  "truncation  error". 

(ii)  In  any  practical  calculating  device,  the  results 
of  arithmetic  operations  must  be  rounded  to  a  finite  number 
of  decimal  or  binary  places.  Thus,  we  are  concerned  also 
with  "round-off  error." 

Obviously,  the  y^  and  y  on  the  right  of  (4.1.3) 
will  not  be  known  exactly,  since  they  themselves  will  In 
general  have  been  calculated  by  (4.1.3)  at  an  earlier 
stage,  and  will  suffer  from  truncation  and  round-off  errors 
The  study  of  the  accumulation  of  errors  assumes 
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crucial  importance  as,  we  shall  see  later.  We  have 
already  had  a  taste  of  the  problem  in  Section  3-5- 


We  begin  the  derivation  of  formulas  by  obtaining 
two  forms  of  Newton’s  interpolation  formula  in  terms  of 
backward  differences.  Regarding  p  as  continuously  variable, 
we  write 


f 


n+p 


(4.1.4) 


Now,  by  definition, 

V  =  1-E 


and  therefore, 

E  =  (1-V)-1, 


(4.1.5) 


whence 


EP  =  (l-V)  p 


(4.1.6) 


=  1+pV+p ( p+1 ) .  p(p+l)(p+2)  V^+. • . 

2 !  31 


(4.1.7) 


Also , 

Ep  =  E(E  P_1 ) ,  (4.1.8) 


and  thus 


EP=E 


1+(P 


1)v+h-l!i£lvg+ 


(4.1.9) 


Thus , 


applying  (4.1.7)  and  (4.1.9)  to  y  n+p. 
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we  obtain 


y'  =  y'+pVy'  +  B(P+1)  v2y'+. 
^n+p  ^  yn  21 


+ 


P  (  )  •  •  ♦  (  P  +  b  1  )  *  _|_E 

kl  ^n  p* 


(4.1.10) 


where 


E  =  hk+1  )  •  ♦  •  (p+^Opk+gy^)  . 

p  (k+l)> 


(4.1.11) 


and  also, 

4+p=  yn+l+<P-1)Vy41+(j4fhi  V2yhn  +  ... 


n+1  21  Jn+1 


+ 


(p _ 1 ) ( P ) ( P+1 ) • « . ( P+k  2 )  yk  1  +E* 

k!  ^n+l  p* 


(4.1.12) 


whe  re 


Ep  =  hk+1  (Pz2g(^+l)|  ;  •  (P^-1)^2^*)  ,  (4.1.13) 


and  where  both  r\  and  i]*  are  in  [xQ,xn+ph]. 

Numerical  integration  formulae  are  obtained  by 
substituting  (4.1.10)  or  (4.1.12)  into  the  relation 

yn+l=  yn-r+h/_r  yn+pdp’  (4.1.14) 

to  obtain,  respectively,  forward  integration  or  iterative 


formulae . 


- 
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For  example,  putting  (4.1.10)  into  (4.1.14) 
with  r  =  0,  we  have 

*^n+l=^n+h(  1+i"  ^+12"  •  •  •  )yn+^n'  (4.1.15) 

which  is  the  Adams  method.  Taking  the  special  case  where 
all  the  differences  are  discarded,  we  have 

yn+l=yn+hyn+h2y"  (4.1.16) 


where  rj  €(xn*xn+]_)*  This  is  yet  another  derivation  of 
Euler’s  method,  this  time  including  an  error  term. 


Rather  remarkably,  if  we  substitute  (4.1.10)  into 
(4.1.14)  and  keep  r3th.  differences  for  r  odd,  the  coefficient 
of  the  r'th.  difference  is  zero  [16],  and  so  r-1  or  r^h  dif¬ 
ferences  provide  the  same  accuracy.  For  r=l,3>5>  we  have, 
with  all  difference  terms  expanded, 

3 

yn+l=yn-l+2hyn+  f  y'"^^  (4.1.17) 


and 


yn+l=yn  -3+311  ( 2yn'yn-l+2yn-2  >  M  ’  ( 4 ' 1 ' 18  > 

y  n=y  ,-+^7r.  (lly  -l4y  n+26y  0-l4y  0 

^n+l  Jn-5  10  v  Jn  Jn-1  Jn-2  Jn-3 


+1 -4 ) +T5oh?yv  1 1  ( 9 )  > 


(4.1.19) 
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where  rj  e  (x  ,x  ,  n  )  . 
1  v  n-r*  n+1' 


Note  that  none  of  these  is  self-starting,  since 
they  all  depend  upon  data  at  least  two  steps  back. 


Now,  putting  (4.1.12)  into  (4.1.14),  with  r  =  0, 
and  keeping  only  one  difference,  we  have 


h. 


3 


w  1  ,  1  \  h.'J  "■/  n 

yn+l=yn+2(yn+l+yn)-  l2y 


(4.1.20) 


with  rj  6  (xn,xn+i  )  * 

Obviously,  we  cannot  employ  (4.1.20)  directly, 

t 

since  we  have  no  value  for  yn+-j_  available.  We  may,  however 
proceed  in  the  following  way: 

(i)  Obtain  some  rough  estimate  yn+1  of  the  value  of 
yn+1 .  For  instance,  (4.1.16)  or  (4.1.17)  might  be  used 
to  obtain  an  estimate. 

(ii)  Using  yn+1  directly  in  the  differential  equation, 

_  t  3 

(4.1.1),  obtain  an  approximate  value  yn+1  of  yn+1 • 

(iii)  Using  (4.1.20),  with  the  estimate  yR+1 ,  obtain 
a  more  accurate  value  of  yR+^ . 

(iv)  Repeat  steps  (ii)  and  (iii)  until  two  successive 
values  of  yn+^  agree  to  the  desired  accuracy. 

Usually,  one  would  choose  h.  small  enough  to 
provide  convergence  after  very  few  iterations. 
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The  use  of  one  approximation  (in  this  case 
(4.1.17) >  called  a  "predictor"  followed  by  Iterated  use 
of  another  (in  this  case  (4.1.20)),  called  a  "corrector" 
is  generally  called  a  "predictor-corrector"  technique. 

The  formulas  (4.1.17)  to  (4.1.18)  are  predictors, 
of  the  second,  fourth,  and  sixth  "order"  respectively. 

The  "order"  of  any  of  these  formulae  is  defined  to  be  one 
less  than  the  power  of  h.  in  the  error  term. 


Corresponding  to  (4.1.18)  and  (4.1.19)  we  may 
obtain  the  corresponding  fourth  and  sixth  order  correctors. 
Again  substituting  (4.1.12)  into  (4.1.14),  it  can  readily 
be  shown  that,  for  odd  r,  the  coefficient  of  the  (r+2)th. 
difference  is  zero,  and  that  we  need  therefore  consider 
only  (r+l)  differences.  For  r=l  and  3*  respectively,  we 
have 


y. 


n+1' 


h  /  1  ,  1  3 

=y  n +^(y - . , +4y  +y 
Jn-1  3VJn+l  Jn  Jn- 


■1 


)-  |o  yV(4) >  (4.1.21) 


and 


yn+l=yn-3+lt(  7yn+l+32yn+12yn- 1 


+32yn-2+'^yn-3) ”555  h7yV11(il), 


(4.1.22) 


where  r\  is  in  (xn_i*xn+i)  and  (xn-3'xn+l^  respectively. 
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Note  that  the  numerical  coefficients  of  the  error 
terms  of  (4.1.21)  and  (4.1.22)  are  considerably  smaller 
than  those  of  (4.1.18)  and  (4.1.19)*  respectively.  This 
is  the  basic  motivation  for  using  a  predictor-corrector 
technique  rather  than  a  predictor  alone. 

An  initial,  obvious,  disadvantage  of  the  above 
forward  integration,  and  predictor-corrector  techniques, 
is  that  they  are  not  self -starting .  That  is,  at  the  be- 

i 

ginning,  to  obtain  y  ,  we  need  not  only  yQ  and  yQ  but  also 

i  t 

y_1,y_0...,  and  y_1*y_2>---  •  In  practice,  we  could  use 

a  Taylor  series,  Euler1 s  method  with,  very  small  h,  or  a 
Runge-Kutta  method  (to  be  described  in  Section  4.4)  to  get 
the  solution  going. 

We  discuss  the  practical  considerations  of  a 
particular  predictor-corrector  procedure,  Milne’s  Method, 
in  the  next  section. 

Section  4.2  Milne’s  Method 

Equations  (4.1.18)  and  (4.1.21)  form  the  basis 
of  Milne’s  method  [ 18 ] .  For  completeness,  we  write 
Predictor : 

—  4h /  ’  *  s  x 

yn+l=yn-3+3(2yn-yn-l+2yn-2)  ’ 


(4.2.1) 


, 


■ 
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and 


_ i 

y 


n+1  1  ^  xni-l  ’  ^n+l  ^ J 


(4.2.2) 


Corrector : 

^n+l_^n-l f 3  ^ yn+l+i+yn+yn-l )  .  (4.2.3) 

The  predictor  and  corrector  are  to  be  used  as 
described  in  the  previous  section;  that  is,  approximate 

_  _ 3 

values  yn+1  and  yn+1  are  obtained  from  (4.2.1)  and  (4.2.2) 
and  then  successively  improved  by  repeated  use  of  (4.2.3) 
and  (4.2.2) . 


The  truncation  error  at  each  step  can  quite 
easily  be  estimated.  According  to  (4.1.18),  the  truncation 
error  of  (4.2.1)  is 

Tl=55‘h5yV(T1l) '  (4.2.4) 

for  e  (xn_3  *xn4  ^)  •  Similarly,  the  truncation  error  of 

(4.2.3)  is 

t2=  -gg-hV^f,  (4.2.5) 

for  n2  Uxn_1<Xn+1). 

Assuming  the  continuity  of  yv(x),  we  may  thus 


write 
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yn+l  -yn+l=95h5yV  ( ^2  >  ^  V  ( ' 'll ) 


(4.2.6) 


for  some  rj  e  (xn_3,xn+i )  *  ff  we  further  assume  that  yv(x) 
is  relatively  constant  in  (xn_^,xn+1) ,  then  approximately. 


yn+l  yn+l" 


29h5  v 
90  y 


(T]): 


29m 

251 


(4.2.7) 


-29T, 


Thus,  by  calculating  yn+1~yn+1  each,  step  we 
obtain  a  good  indication  of  the  local  truncation  error.  If 
this  quantity  grows  smaller  over  several  successive  steps, 
then  we  may  increase  h.,  or,  if  it  grows  larger,  we  must 
decrease  h. 

One  major  difficulty  with  the  method  described 
above  is  that,  in  the  process  of  iterating,  we  must  evaluate 
the  function  f  each  time  around.  Thinking  again  of  sets 
of  equations,  we  realize  that  f  will  often  be  quite  complicated 
to  evaluate,  and  that  a  good  deal  of  computing  time  may  thus 


be  involved  in  each  iteration. 


- 
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A  modification  of  Milne’s  method  which  counters 
this  objection  by  requiring  only  two  evaluations  of  f  at 
each  step  can  be  used,  at  a  modest  cost  in  accuracy. 

We  change  our  notation  slightly  and  write 

Pn+V  °n=V  (4-2'8) 

Thus,  we  rewrite  (4.2.1)  as 

Pn+l“^n-3+5^24i~'4-l+25rn-2)  ’  (^-2.9) 


Also  (4.2.7)  becomes 


cn+l  ^n+1" 


2  SL  5  v 


90 


■h  y 


(4.2.10) 


Let  us  assume  that  yV(r])  is,  to  a  fair  approximation, 
constant  from  one  step  to  the  next.  Then,  using  values 
from  the  n’th  step,  we  would  find,  from  (4.2.10), 


m  28  /  \ 

Tn  ~  c  -p  )  , 
1  29  n  ny 


(4.2.11) 


and 


so,  adding  T^  to  Pn+1  should  give  a  very  good  estimate 


of  c  , ,  .  We  write 
n+1 


pQ 

m  .  =  p  ,n  +  7Trr( c  “P  )  ) 

n+1  n+1  29  n  n7 


(4.2.12) 


, 
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and 


"n+l  =  f(xn+l'ran+l)- 


Now,  using  mn+1,  we  rewrite  (4.2.3)  as 


°n+l  =  yn-l+J(mn+l+4yn+yn-l} 


(4.2.13) 


We  can  further  improve  this;  again,  from  (4.2.7),  we  may 
write 

T2  %  29  ^Cn+l“Pn+l^,  (4.2.14) 

and  so,  adding  to  c  we  have,  finally. 


yn+l  cn+l+  29^Pn+l  cn+l^ 


(4.2.15) 


and 


yn+l  =  f^n+l’^n+f- 


(4.2.16) 


The  quantities  p  are  called  the  "predictors",  m 
the  "modifiers"  and  c  the  "correctors." 

This  scheme  involves  the  evaluation  of  f  only 
twice  at  each  step,  and  is,  in  practice,  very  nearly  as 
accurate  as  the  iterative  scheme  using  (4.2.1),  (4.2.2), 
and  (4.2.3)-  Furthermore,  we  can  use  the  quantity 


-  f,,fr 


..  .. 


>9  J'; 


J  !  .  ‘  '  £r .;  ••  [V: 


999  ".i ,  l  n  ,o  oVl 

9  j  :'i  / 

i  .  v<- 


,v.;.".oa.rr  ,  9V£(1 


'  OB;  O/IB 


s  l~ 


[  r+i  "'**!'■  q  ’  >r  '  .  4-n 


1  .  "  -  n  '  .  Li  L-  ■  ■■  •  .1  ■  i'  .'i  I...'  9  i  ; 

1  ..  .■  "o  ■  oo '.to:  ..  : ;  9:L;  o  r :o\  :l :  •:  ;  ,:t' doav"  .  L; 


■  ..  ■  ;  1 .  9 "vv  !'  •  '■  .  i  ,9Lr .  99  ,  9':.' 


,;Im  o:.'  ■;  3SV  t -V  LlLL.CL  ffl  ,  .  !:';i..!.k  .  i  ftiOSQ  JL  30.;  Wt 

-V 

(  ,  c  .  (  J  .  )  :■  ,9;  319..  L  .  f  -  Ct,.  ”  •  It  c  ■  9  £'  9 

J93  if. ' :j  ,  9  9"  *  ,  ‘9  , ;  '  .  '  .  !  g  ..  ■' ,  - 


D11B 


90 


c  +i"Pn+i  as  an  estlmate  of  the  error  at  each  step. 

At  this  point,  we  return  to  the  original  problem 
of  this  section  and  describe  this  process  algorithmically 
for  the  set  of  first  order  equations  ( 3 •  3 •  )  subject  to 
(3.3*6)  .  To  do  so,  we  represent  the-  vectors  ^-n-2’ 

and  as  the  columns  Y-  ,  Y^,  Y^,  and  Y^  respectively 

i  i  i 

of  a  matrix  Y.  The  vectors  y  y  , ,  and  y  we  represent 
as  the  columns  W^,  W^,  and  of  a  matrix  W.  Also,  the 
vectors  p^  and  are  represented  as  the  columns  P-^  and 

P_2  the  matrix  P.  ?  is  a  zero  vector.  It  is  assumed 
that  some  starting  process  has  calculated  and  output  the 
values  of  Y^,  Y^,  Y^,  Y^,  W^,  W^,  and  W^,  for  x  =  xq, 
xQ+h,  xQ+2h.,  and  xQ+3h.  respectively.  The  vectors  m,n  and 

i 

c_  correspond  to  the  quantities  mn+^,  and  cn+l 

equations  (4.2.12),  (4.2.12),  and  (4.2.13)*  respectively.: 

MILNE  1 S  METHOD 


1 

2 

3 

4 

5 

6 


x  «—  4h.  +x 


o 


4h./ 


initialize  x 
initialize  c 


initialize  _P  ^ 


P_2  «—  Y -l+~^(  2W2-W2+2W-l  )  calculate  predictors 


—  4  -2+29^-  -1^ 
n  «—  f (x,m) 


calculate  modifiers 

calculate  derivatives  of 
modifiers 


0<O 


i+fl  '  i  .0 


1 1 :  :  o '  ;  oT 


Qt  do  'c 
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7 

c  4-  Y3+~(n+4W3+W2) 

calculate  correctors 

8 

Shift  columns 

9 

-2  -3 

y  of  Y  left 

10 

J 

one  position. 

ll 

^  Wg 

1  Shift  columns  of  W 

12 

W2  ^  w3  , 

left  one  position. 

13 

X4  — +29^— 2“— ^ 

calculate  final  values 

14 

«3  - 

calculate  final  derivatives 

15 

Pn  P„ 

— 1  — 2 

shift  columns  of  P  left. 

16 

PRINT:  x, 

output  results 

17 

x  «—  x+h 

increment  h. 

18 

x:x  *(<*>)— *-(4, 19) 

test  for  final  x 

19 

STOP. 

Inspection  of  this  program  shows  that  we  are 
carrying  four  columns  of  Y,  three  of  W,  two  of  P_,  and  one 
each,  of  c_,  m,  and  n,  for  a  total  of  12 1  computer  storage 
locations,  where  H  is  the  number  of  simultaneous  equations 
being  solved.  Acutally  steps  5*6,7  could  be  written,  some¬ 
what  awkwardly,  as  one,  thus  eliminating  the  vectors  m  and 
n.  Hence,  we  require  10i5  storage  locations. 

It  is  possible,  using  c  -  as  a  criterion  to 
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gauge  the  truncation  error  at  each  step.  Based  on  this, 
we  can  increase  or  decrease  the  interval  h.  to  gain  speed 
when  possible,  or  to  preserve  accuracy  when  necessary. 

Such  a  scheme  is  described  in  Ralston  and  Wilf  [16],  but 
requires  storage  for  additional  columns  of  Y  and  W. 

The  routine  described  above  should  be  fairly 
fast,  the  major  time  being  spent  in  the  two  evaluations  of 

f . 


Thus,  Milne’s  method  seems  ideal,  except  for 
one  thing.  Nowhere  in  Milne’s  original  paper  [18],  and, 
indeed,  nowhere  in  his  later  book  [20] ,  is  there  any  dis¬ 
cussion  at  all  of  stability.  In  fact,  the  word  does  not 
appear  to  be  used  at  all  in  either  reference. 

As  we  shall  see  in  the  following  section,  Milne’s 
method  suffers  from  definite  instability,  but,  it  can  be 
modified  at  fairly  small  cost,  to  guarantee  stability  under 
fairly  mild  conditions. 

Section  4.3  Stability  of  Multi-Step  Methods 

Recalling  the  example  of  Section  3.5>  we  saw  that 
a  component  of  an  undesired  solution  of  a  difference 
equation  could  grow  and  dominate  the  desired  solution.  In 
that  section,  we  were  approximating  a  second-order 
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differential  equation  by  a  second-order  difference  equation. 
With  the  multi-step  methods  of  the  previous  two  sections, 
we  have  been  approximating  first  order  differential  equations 
by  higher  order  difference  equations.  Accordingly,  we  must 
be  extremely  wary  lest  one  of  the  several  solutions  of  such 
higher  order  difference  equations  misbehave  and  dominate 
the  desired  solution. 

In  this  section,  we  present  an  analysis  of  the 
stability  of  Milne’s  method  by  a  method  which  is  of  general 
applicability. 

We  will  consider  the  iterative  version  of  Milne’s 
method,  since  the  analysis  is  somewhat  simpler  than  the 
non-iterative  method,  and  the  stability  properties  of  the 
two  are  not  greatly  different. 

Following  Hamming  [21],  we  consider  the  "general¬ 
ized  corrector"  formula 


yn+l  =  Pyn  +  qyn-l  +  ryn-2 


+h(syn+l  +  ty^  +  uy^),  (4.3.1) 


in  place  of  (4.2.3). 
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Let  z  be  the  true  solution  of  the  given  problem 


that  is 


1  dz  „/  \ 

z  =  57  =  f(x’z)’ 


(4.3.2) 


whereas  the  solution  which  we  calculate  satisfies 


y  =  f (x  ,y  )  +  E(n) , 
Jn  v  n  Jn'  '  '  ’ 


(4.3.3) 


where  E(n)  is  the  error  at  the  n'th  step,  and  is  supposed 
to  be  small. 

Thus,  we  can  write 


z  ,,  =  pz  +  qz_  ,  +  rz  0 
n+1  ^  n  ^  n-1  n-2 


+h(  szn+l.  +  tzn  +  uzn-l)+F(n)'  (^.3. 


where  F(n)  is  the  corresponding  error. 


We  define  £,  the  error  of  the  solution,  as 


e  =  z  -  y  . 
n  n  Jn 


(4.3.5) 


Subtracting  (4.3.1)  from  (4.3.4),  we  have 


£  , ,  =  p£  +  q£  +  r£  ~ 
n+1  ^  n  H  n-1  n-2 

+h^sen+l  +  t£n  +  +  F(n) 


'n-1' 


(4.3.6) 
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Also,  subtracting  (4.3.3)  from  (4.3.2),  we  have 

en  =  f(xn,zn^"f^xn,yn^"E^n^ ‘  (4.3.7) 

Now,  applying  the  mean-value  theorem  to  (4.3*7) > 

en  =  Sf(Xn,Tii7  e  -E(n),  (4.3.8) 

n  Zy  n 

for  some  n  e [y  , z  ] . 

'n  °n  n 

In  practice,  E(n),  F(n),  and  ^  change  fairly 
slowly  from  step  to  step,  and  so  we  will  treat  them  as 
constants . 


Writing 


A  = 


df 
dy  ' 


(4.3.9) 


we  transform  the  independent  variable  x,  to 


x  =  t/A. 


(4.3.10) 


In  these  units,  we  may  now  rewrite  (4.3.8)  as 


d£  _  __  E 
dt  “  A’ 


(4.3.11) 


which  we  substitute  into  (4.3.6),  to  obtain 


£  . ,  =  p£  +  q£  ,  +  r£  0 
n+1  n  ^  n-1  n-2 


+  k( s£  ,n  +  ts  +  u£  n ) 
v  n+1  n  n-17 


+  F  -  ^-(s+t+u)E, 


(4.3.12) 


where  k  =  hA. 

Equation  (4.3.12)  is  a  third-order  linear  difference 
equation,  which  we  attempt  to  solve  using  the  trial  solution 


e  =  p 
n  r 


n 


(4.3.13) 


This  yields  the  indicial  equation 


=  Pp2+qp+r+k(sp3+tp2+up) ,  (4.3,14) 


or 


3  2 

k  =  p  -Pp  -qp-r 
3  2 

sp-^+tp  +up 


(4.3.15) 


For  Milne’s  method,  we  evidently  have 


1  4 

p=r=0,  q=l,  s=u=j,  t=j. 


(4.3.16) 


' 


■ 


whence  (4.3.15)  becomes 


(4.3.17) 


which,  for  any  given  k  has  two  roots  p^  and  p^. 

Thus,  the  general  solution  of  (4.3.11)  is 


e 


n  .  n  , 

n  ~  clpl  +  2P2  +  7 


(4.3.16) 


for  some  y. 


Hamming  shows,  by  a  graphical  technique,  that  one 


or  the  other  of  p^  of  p^  is  always  greater  than  1  in 
absolute  value  whether  A  is  positive  or  negative.  Hence 
Milne 3 s  method  is  definitely  unstable. 


We  point  out  also  that,  in  practice,  the  absolute 


magnitude  of  the  error  is  not  the  prime  consideration; 

rather,  we  are  usually  more  interested  in  the  size  of  the 

error  relative  to  the  size  of  the  solution.  Thus,  we  are 

led  to  the  concept  of  "relative  stability,"  implying  that 

the  relative  error  of  an  approximation  remains  bounded. 

For  relative  stability,  we  may  be  more  interested  in  knowing 

_h 

that  some  function  such,  as  pe  '  is  less  than  unity  in 
absolute  value. 
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In  addition  to  demonstrating  the  instability 


of  Milne's  method,  Hamming  shows  an  acceptable  alternative 
procedure.  He  observes  that  Milne's  scheme  is  exact  for 
a  polynomial  of  degree  four  or  less. 


Hamming  imposes  the  auxiliary  conditions  that 

o  o 

the  corrector  formula  (4.3.1)  be  exact  for  y=l,x,x  ,x  ,x  . 
Thus,  he  obtains  the  relations 


p  =  27(l-q) 
24 


q  =  q 


s  =  2-za. 
s  24 


r  = 


-3(i-q) 

24 


t  = 


u  - 


l8+l4q 

24~^ 


-9+17  q 

24 


(4.3.17) 


j 


He  further  shows  that  the  error  term  is 


mh5yV(r])  =  h5yV(r] )  .  (4.3-18) 

For  h.  =  k  =  0,  (4.3*13)  becomes 

p3-pp2-qp-r  =  0,  (4.3.19) 

or,  using  the  values  from  (4.3.17)* 

8p3-9(l-q)p2-8qp+(l-q)  =  0.  (4.3-20) 
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Plotting  (4.3.20)  as  a  graph  of  p  versus  q, 
Hamming  shows  that  the  value  q=0  (for  Milne's  method, 
q=l)  provides  a  good  compromise  between  the  requirements 
of  stability,  and  of  ease  of  evaluation. 

The  method  implied  by  q  =  0  is  often  called 
"Hamming's  Method."  It  is  stable  and  relatively  stable 
for  A<0  if 

h<  ,  (4.3.20) 

lAl 

and  is  relatively  stable  for  A>0  if 

h<  —  *  (4.3.21) 


A  variation  on  Hamming's  method,  in  the  same 
vein  as  the  version  of  Milne's  method  which,  does  not  require 
iteration,  requires 

h<  & 

A 

This  latter  method  may  be  represented  algorithmically 
by  replacing  steps  5 >7 ,  and  13  in  the  algorithm  for  Milne's 
method  by 
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112 

121 
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5-(9Y4-Y2)+3h.(n+2W3-W2) 
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121 
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The  procedure  described  above  for  analyzing 
the  stability  of  the  Milne  process  may  be  employed,  with, 
obvious  modifications,  to  the  analysis  of  many  multi-step 
procedures.  Detailed,  and  more  general  analyses  are  given 
by  Urabe  [22],  and  by  Hull  and  Luxemburg  [23]. 

Section  4.4  Runge-Kutta  Methods  for  Initial  Value  Problems 

The  methods  of  the  previous  three  sections 
are  called  multi-step  mthods  because  they  make  use  of 
information  several  steps  back  in  the  integration  process. 
In  this  section,  we  consider  several  members  of  the  class 
of  one-step  methods,  so-called  because  they  make  no  use  of 
previous  data. 

We  recall  that  one  of  the  major  objections 
to  the  Taylor  series  method  was  that  expressions  for 
derivatives  of  f  are  difficult  to  calculate.  The  Runge- 
Kutta  methods  avoid  this  difficulty  by  using  additional 
evaluations  of  f  itself,  to  attain  comparable  accuracy. 

Following  Henrici  [9] >  we  demonstrate  the 
basic  approach  by  developing  a  Runge-Kutta  formula  of  the 
second  order  for  the  usual  problem. 

y  (x)  =  f(x,y) , 


(4.1.1) 
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subject  to 


y(x0)  =  y0 


(4.1.2) 


By  the  Taylor  series  method,  we  can  write 


2  m 


y(x+h)-y(x)  =  hy* (x)  +  |y  y  (x)  +  0(h3),  (4.4.1) 


or. 


dx+hJbylx),  =  y  (x,y;h). 


(4.4.2) 


where 


T ( x , y j h. )  -  f(x,y)+|f '  (x,y)+0(h2)  (4.4.2) 


By  way  of  an  approximation  to  T(x,y;h),  let 


us  consider 


<t>(x,y;b)=  ay  (x,y)+a2f  (x+p  h,y+p2hf(x,y) ) , 


(4.4.3) 


with,  the  parameters  a^a^^p^  and  p^  at  our  disposal. 
Expanding  in  a  Taylor  series,  we  have 


f  (x+p1h.,y+p2hf  (x,y) ) 


f(x,y)+P1h.-f^y^+p2h.f(x,y)^-^^+0(h2) , 


(4.4.4) 


■ 
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and  thus  we  may  rewrite  0  as 


$(x,y;h)  =  (a1+a2)f(x,y)+a2h. 


df(x,,y) 


1  dx 


+  P2f(x,y)^|^  +0(h2).  (4.4.5) 


Also,  by  elementary  calculus,  and  using  (4.4.1), 


d 


f  =  dxf^x,y^ 


Sf^x)+  f(x,y) 


(4.4.6) 


and  thus  we  may  rewrite  'f  as 


S(x,y;h)  =  f(x,y)+' 


h 


^y)+f(x,y)^z) 


+0 ( h2 ) 


(4.4.7) 

2 

We  now  require  that  0  and  Y  agree  within  0(h  ) . 
Thus,  matching  terms  in  (4.4.5)  with  (4.4.7)*  we  have,  from 
the  constant  terms. 


al+a2  =  ^  ’ 


(4.4.8) 


and  from  the  coefficient  of  h., 


a2pi  2’  a2P2  2  * 


(4.4.9) 


•  I  J  Mi  ,  •  < 
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Thus,  writing  a  =  an  we  have  immediately, 

a1  =  1  -cl,  a2  =  a,  p1  =  P2  =  ^  *  (4.4.10) 

and  a  is  free  to  take  on  any  non-zero  value.  Hence,  we 
write  0  as 

$(x,y;h)  =  (l-a)f(x,y)+af {x+^y+^f (x,y) },  (4.4.11) 

2 

and  note  that  0  deviates  from  T  only  by  0(h.  )  for  any 
value  of  a. 


Two  particular  values  of  a  are  historically 

important : 

(i)  For  a  =  l/2,  we  have  the  "improved  Euler  method" 
or  the  "Heun  method."  The  approximation  is  written 


y  ,  =  y  +^-{f(x  ,y  )+f(x  +h,y  +h.f(x  ,  y  ))); 

°n+l  n  2  n  ^n'  n  Jn  v 

(4.4.12) 

(ii)  For  a  =  1,  we  have  the  "improved  polygon  method," 
or  "modified  Euler  method,"  written 

h.  ,  h. 


y  ,i-y  +hf(x  +7T, y  +7rf ( x  ,y  )). 
Jn+1  ^n  vn2n2nJn// 


(4.4.13) 


Notice  once  again  that  both  of  these  formulas 


have  the  same  order  of  truncation  error  as  the  Taylor  series 


' 
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method  with  one  derivative  term,  but  that  they  both 
achieve  this  by  two  evaluations  of  f,  and  no  evaluations 
of  the  derivative  of  f. 

The  classical  Runge-Kutta  formula  [24,25]  is 
similar  to  the  one  we  have  discussed,  but  it  employs  a 
weighted  sum  of  four  evaluations  of  f,  and  thus  has  the 
same  order  of  truncation  error  as  a  Taylor  series  with,  the 
third  derivative  term  of  f. 

Specifically, 

^(x,y.;h)  —  (-L]_k2.+^2k2+tl3k3+^4k4  (4  4  14) 


where 


k1=  f (x,y) 

k2  =  f  (x+ah^y+pk^) 

k^  +  f  (x+a1h.,y+51k1+71k2) 

k^  =  f (x+a2h,y+p2k1+72k2+52k^) 


(4.4.15) 


Relations  among  the  aJs,  5’s,  yJs,  &2,  and  \± 1  s 
are  obtained  by  requiring  the  Taylor  expansion  of  (4.4.14) 
to  agree  with  the  Taylor  series  method  up  to  terms  of 


' 

■ 


■ 
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order  h  .  The  algebra  Involved  in  this  process  Is  quite 
lengthy,  but  Involves  only  standard  techniques,  and  is  of 
no  particular  interest  in  itself;  accordingly,  we  omit  it 
from  this  discussion.  The  details  may  be  found,  for  example, 
in  Ralston  and  Wilf  [ 16 ] . 


This  process  leads  to  eight  equations  in  the 
thirteen  parameters.  By  requiring  these  to  be  independent 
of  f  itself,  one  obtains 


a  =  (3 

al  =  ^1+71 
a2  =  ^2+72+^*2 


\ 

► 

j 


(4.4.16) 


and  the  following  eight  equations  in  ten  parameters 


*J'1+IJL2+  ^3+^4 

2  2  2 

\^2a 

+q^a2 

+b^(a72  H-a^) 

2  2  2 
l-^a  7]_  +M-4(a  72+a162^ 

M’3aa171+,J'4  ^ a72+al62  ^ a2 

^4a71^2 


=  1 

=  1 
2 

=  1_ 
3 
1 

=  4 


=  1 

=  _JL 
12 
=  1 
E 
1 

=  24 


> 


(4.4.17) 


, 


. 


'  ■  , v  '  ■"  :  1  .  :  V.  : 
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Equations  (4.4.17)  have  two  degrees  of  freedom; 
most  of  the  variations  of  the  Runge-Kutta  scheme  make  use 
of  this  freedom  in  slightly  different  ways. 


Runge  [24],  used  the  following: 


=  ^4  =  ) 

m«2  =  ^3  ~ 

kx  =  f (x,y) ,  : 

k2  =  f ( x+h/2 , y+kx/2 ) , 
k^  =  f  (x+h./2,y+k2/2) , 
k^  =  f(x+h,y+h).  j 


Kutta  [25]  *  chose  instead: 


M-i  =  M-4  =  V8* 

M-2  =  M3  =  3/8, 

k-L  =  f(x,y), 

k2  =  f  ( x+h/3 ,  y+^/3 ) , 

k^  =  f  (x+^,y-k1/5+k2) , 

k^  =  f  (x+hjy+^-kg+k^)  • 


(4.4.18) 


(4.4.19) 


These  choices  appear  to  have  been  motivated  mostly 
by  a  desire  to  make  the  algebra  as  simple  as  possible.  More 
recently.  Gill  [26] ,  motivated  by  the  availability  of  a 
digital  computer  (EDSAC)  with,  a  small  memory,  and  also  by 


, 
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a  desire  to  minimize  the  round-off  error,  used  the  following, 
which  is  commonly  known  as  the  Runge-Kutta-Glll  procedure, 
and  is  probably  the  most-used  method  for  solving  initial  value 
problems  using  a  computer: 


kl  =  f(xo’yo}’ 
k2  =  f(xQ+h/2,y1), 

k3  =  f(x0+h/2>y2) ' 
k4  =  f(xo+h,y3); 


yl  =  yo+(kl-2qo)/2' 

y2  =  y1  +  (lWl/2)  (k2-q1)  , 

y3  =  y2+(iWi/2)  (k3-q2) , 

y^  =  y3+(kZ|-2q3)/6; 


ql  qo+2^kl  2qo>  kl/2* 

q2  =qx  +3(1-^172)  (k2-q1)-(lVl72)k2^ 

q3  -  q2+3(lWI72)(k3-q2)-(l-iVl72)k3, 

O  1 

G4  =  q3+§(k4-2q3)-2k4; 


7  (4.4.20) 


y 


where  y(xQ+h)  =  y^,  and  qQ  =  0  initially,  and  qQ  =  q^ 
after  the  first  step. 


Clearly,  as  each  of  the  k*s,  y's,  and  q5s  is 
calculated,  it  may  replace  its  previous  value. 


We  represent  the  Runge-Kutta-Gill  process 
algorithmically  for  the  set  of  equations 

y' (x)  =  f(x,jr) , 


(3.3.4) 


9  ■  ..  :>•  ■  ■'  .  ,  'V  -  •  $  fj:;  ;  '  v  Licon  :•  .  ;)i  it)  UU 
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subject  to 

z(xQ)  =  a-  (3.3.6) 

The  vectors  y,  and  a  represent  the  obvious 
quantities  of  (4.4.20),  and  $,  Is  the  number  of  equations. 
We  assume  that  the  auxiliary  vectors  a,  b,  and  c_  are 
Initially 

a  =  (l/2,  1-41/2,  1+45/2,  1/6), 
b  =  (2, 1,1, 2), 

c  =  (1/2,  1W172,  1WI72,  1/2). 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 


RUNGE-KUTTA-GILL  METHOD 


x  «—  x 

o 

t  <-  a 

£  «-  S_ 

j  -  1 

k  <-  f  (x,y) 
i  «—  1 

z  «—  a  . (k.  -b  .q.  ) 

yi”yi+hz 

q.x-q .  +  3  -c  .k. 

^1  z  j  1 

i  •<—  i+1 

j*-j  +  l 


initialize  x 

initialize  y 

initialize  £ 

initialize  j 

evaluate  k 

initialize  i 

calculate  auxiliary  z 

calculate  new  y. 

calculate  new  q. 

^•1 

increment  i 
test 

increment  j 


,  ■■  .  .  0  J  ’T 


;  v  i  ■  c 
i. 
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13  j:4,(£,>M6,l4) 

14  PRINT;x,y 


output  solution  point 


test 


15  x«— x+h 


Increment  x 


16  x:x 


'max 


(i>>M4,17) 


test 


17  STOP 


We  notice  several  things  about  the  Runge-Kutta- 


Gill  process,  by  way  of  comparison  with  Hamming's  method. 
One  great  advantage  of  the  Runge-Kutta  scheme  is  that  it 
is  self -starting,  unlike  the  Hamming  procedure.  On  the 
other  hand,  it  requires  four  evaluations  of  f  for  each 
step  forward,  and  will  therefore  require  about  twice  as 
long  as  Hammings  method,  on  a  computer.  Although  one  of 
Gill's  innovations  was  to  reduce  the  storage  required  to 
3f  locations  (as  opposed  to  10 l  for  Hamming),  this  is  no 
great  advantage  on  modern  computers  which,  usually  have 
large  quantities  of  rapid,  random-access  memory. 

Stability  of  Runge-Kutta  methods  is  seldom 
discussed.  However  the  method  is  sometimes  unstable  (see 
for  example  [11],  page  91] *  although  this  can  usually  be 
overcome  by  taking  h.  sufficiently  small.  A  recent  paper 
by  Urabe  [22]  contains  the  most  complete  treatment  known 


to  us  . 
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It  is  our  experience  that  the  Runge-Kutta 


method  is  over-used  in  many  computing  centres,  on  account 
of  its  ease  of  use.  We  would  suggest  that  a  much  better 
general-purpose  program  would  merely  use  a  Runge-Kutta 
starter  for  a  Hamming  process,  for  the  reasons  which  we 
have  discussed  above. 

Section  4.5  Further  Techniques  for  Boundary  Value  Problems 

We  again  consider  the  linear,  second  order 

equation 


y  (x)+f (x)y' (x)+g(x)y(x)  -  k(x),  (3.6.3) 


ii 


subject  to 


y(a)  =  a,y(b)  =  p. 


(3.6.2) 


Using  the  approximations 


2  _4  c6 


(2.3.11) 


(2.3.15) 


(3.6.3)  becomes 


(4.5.1) 


. 

'  0  ■  .  1:  '  l  l  +  ?.T 


1  1 

Ill 


That  is 

62y1+hfiix5yi+h"giyi  =  h2k1+  Cyi,  (4.5.2) 

where 

c  =  +-.-)+hf1(y:|^  -  ^§4...)  (4.5.3) 

Rewriting  (4.5.2)  in  terms  of  ordinates, 

yi-l(lTfi)"yi(2_h2gi)+yi+l(l+lfi)  =  h2k1+Cy1.  (4.5.4) 


Thus,  in  the  same  manner  as  in  Section  3.6,  we 
can  write  (4.5.4)  as 


A  1=  d-C  y. 


(4.5.5) 


and  the  boundary  conditions  are  included  in  d. 


Now,  for  a  first  approximation,  we  neglect  C_  y, 
and  obtain  a  first  approximation  y^  to  y  by  solving 


A  =  d 


(4.5.6) 


We  next  difference  y^  to  obtain  an  initial 
estimate  of  C  y,  and  calculate  a  correction,  to  y^  by 


*  1  =  -£  jr1. 


(4.5.7) 


' 


.  *0+  jfrt 
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if  c  a1 

further  correction 


is  still  large,  we  can  calculate  a 
,  from 


A 


(4.5.7) 


and  so  on  until  the  corrections  have  become  as  small  as 
we  desire.  Then 

z=  z1+a1+ir+-  •  •  ♦  (4.5.8) 


The  only  difficulty  with,  this  process  is  that, 
in  order  to  evaluate  the  difference  corrections,  using 
central  differences,  we  require  values  of  y  outside  the 
range  [a,b] .  These  can  readily  be  obtained  by  using 
equation  (4.5-4)  as  a  predictor  formula. 

An  example  of  this  process  is  given  in  [11],  page 

95- 


Another  way  to  solve  the  given  boundary  value 
problem  is  to  replace  it  with,  a  pair  of  initial  value 
problems,  as  follows: 


Using  any  of  the  methods  which,  we  have  discussed 
for  solving  initial  value  problems,  obtain  a  solution  y. , 
of  (3.6.3)  subject  to 

y1(a)  =  a,  y[ (a)  =  ■/1, 


(4.5.9) 


- 
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and  obtain 


y1(b)  =  p1 


(4.5.10) 


is 


Then,  solve  the  homogeneous  form  of  (3.6.3)*  that 


y^(x)+f(x)y2(x)+g(x)y2(x)  =  0,  (4.5.11) 


subject  to 


y2(a)  =  0,  y2(a)  =  y2 


(4.5.12) 


Since  the  differential  equation  (3.6.3)  is  linear, 
any  linear  combination  of  y^  and  y^  also  satisfies  it.  In 
particular. 


y  =  yx(x)+ty2(x) , 


(4.5.13) 


not  only  satisfies  (3*6.3)*  but  also  y(a)  =  a,  for  any  t. 


Clearly, 


y(b)  -  y1(b)+ty2(b) 


=  P1+tp2=p* 


(4.5.14) 


and  thus,  solving  (4.5.14)  for  t, 

t  =  P-Pi  . 

6o 


(4.5.15) 


Hence,  the  desired  solution  is 
y(x)  =  yx(x)+  y2(x). 


(4.5.16) 
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CHAPTER  V 


A  NUMERICAL  APPROACH  TO  THE  GREEN’S  FUNCTION 


In  this  chapter  we  consider  the  matrix  approxi¬ 
mation  to  the  boundary  value  problem  given  in  Chapter  III. 

We  show  that,  as  the  size  of  the  matrix  increases,  and  the 
tabular  value  decreases,  the  approximation  takes  on  the 
form  of  the  solution  given  by  the  well-known  Green's  function 
method . 


A  proof  of  convergence  of  the  approximation  is 
given,  and  an  analytic  expression  for  the  Green's  function 
is  obtained.  This  expression  is  believed  to  be  new. 


Section  5>1  Introduction 

In  this  chapter,  we  will  consider  the  linear 
boundary  value  problem 

y"(x)+q(x)y(x)  =  k(x),  (5*1,1) 


subject  to 


y(0)  ^  a,y(l)  ^  p. 


(5.1.2) 


The  differential  equation  (5 • 1 • 1 )  1®  hot  as 
restricted  as  it  might  at  first  appear,  for  we  can  usually 
write 

yn(x)+p(x)y'(x)+q(x)y(x)  a-  r(x),  (5-1*3) 

in  the  same  form  as  (5.1.1),  as  we  now  show.  In  (5.1.3), 
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let 


y(x)  =  u(x)v(x), 


whence 

y  (x)  =  u1 (x)v(x)+u(x)vJ (x) , 


and 

y" (x)  =  u"(x)v(x)+2u'(x)v'(x)+v"(x). 
Thus,  (5*1*3)  may  be  written 

?  J  J 

u'!v+u  (2v  +pv)+u(pv  +qv)  =  r-v  , 
and,  to  eliminate  the  first  derivative  term  in  u, 

i 

2v  +pv  =  0, 

and  we  can  solve  for  v  by  separation  of  variables 
write  (5.1.8)  as 


dv 

dx 


or 


dv 

v 


(5.1.4} 


(5.1.5) 


(5.1.6) 


(5.1.7) 
we  set 

(5.1.8) 

We  c  an 

(5.1.9) 


(5.1.10) 
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and  hence,  integrating  over  a  suitable  range 

In  v  =  -  |  / pdx,  (5.1.10) 

or 

v  =  exp(-  /pdx).  (5.1.11) 

In  this  way,  subject  to  fairly  mild  restrictions, 
equation  (5.1.3)  may  be  written  in  the  form  (5.1.1). 

Let  us  now  solve  (5.1.1)  numerically  by  the 
method  of  Section  3.6.  Obviously,  we  again  have  to  solve 


A  y  =  d , 


(3.6.14) 


where  now  A  is  the  matrix 


(2-h2g1)  -l 

0 

0.  .  . 

-l  (2-h2g2) 

-1 

0  •  .  . 

0  -1 

(2-h2S3) 

-1.  .  . 

•  *  •  t  t  “1 

(2-h2gn_g) 

-1 

V. 

-1 

(2-h26n.i) 

S 


(5.1.12) 


t 


■  ■  -  i  '  a  ,  94- 1 


■  '■  ■  -  1  -  "  '  l-  r  :  .  :,v-  c,  :  i  csi 

c 


.S;*  -  .  .  ;  -  '-f;  j't:  • 


\ 


!>-  '5’ 


J'  C 


*  J?  *  * 
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y  is,  as  before 


and 


d  = 


y- 


y 


n-1 


s  \ 

/  N 

h2k^-a 

kl 

a 

h2k2 

k2 

0 

h2k3 

k3 

0 

• 

CVJ 

A 

II 

« 

" 

• 

• 

• 

• 

• 

* 

t 

h2k  0 

k  ~ 

0 

n-2 

n-2 

h2kn-i  -  e 

kn-l 

P 

k  / 

(5.1.13) 


h2k-b. 


(5-1.14) 


Thus,  the  solution  is  given  by 
y;  =  A”1(h2k)  -A'1^. 


(5.1.15) 


Writing 


(5.1.16) 


. 


■ 


- 


■ 


r 


■ 


, 
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we  have 


y;  =  G(hk)  -  G  b/h. 


or,  for  the  components  of  ^ 


n-1 


n-1 


yi  =  L  Gukjh  '  L  GuVh 

j=l  J=1 


But 


bj  '  a5l,  j+(36n-l,  j’ 


(5.1.17) 


(5.1.18) 


(5.1.19/ 


where 


is  the  Kroneker  delta  function, 


and  so 


n-1 


V  G.  k  .h 

L  u  j 

j=i 


G 

{a- 


1,1 

h. 


+  p- 


G, 

i,n-l 

h. 


)• 


(5.1.2C) 


Assuming  for  the  moment  that  a  =  p  =  0,  we  have 


n-1 

y .  =  )  G.  .k.h. 

1  L  j 

J-l 


(5.1.21) 


Now,  it  is  well  known  [28]  that  the  solution  for  equation 
(5.1.1)  sub j  ect  to  (5.1.2)  with  a  =  p  =  0  can  he  written 


y(x)  =  f  G(x,0lc(0^i 
GO 


(5.1.22) 


’  v  . 


••  7  •  ■  ■  M  -  •  -  .  ■  '  .  i:  ,  ,0  r>r  .y 


■ 


:.v  ■  j  r  -  twtr: 
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where  g(x,£)  is  called  the  Green's  function.  We  point 
out  the  obvious  similarity  between  (5.1.21)  and  (5.1.22); 
although  the  one  is  a  discrete  formula,  and  the  other 
continuous,  it  is  clear  that  the  sum  term  in  the  former 
corresponds  to  the  integral  term  in  the  latter.  In  fact, 
(5.1.21)  looks  like  a  trapezoidal  approximation  to  (5.1.22), 
and  thus,  we  are  led  to  suspect  that,  in  the  limit  as 
h -*■  0,  (5.1.21)  may  coincide  with  (5.1.22);  we  will  show 
that  this  is  the  case  later. 

We  can  go  farther  and  use  (5.1.20)  to  obtain  a 
plausible  extension  to  (5.1.22).  We  remove  the  requirement 
that  a  =  p  =  0;  now  using  the  (so-far  unproved)  correspon¬ 
dence 

G  - >  G(x,5),  (5.1.23) 

1J  h-*0 

we  obtain,  by  taking  the  limit  of  (5.1.21), 

y(x)  -  ffl(x,eW0dHa  “3  +  G(*u-h)). 

(5.1.24) 

Now,  we  know  independently  that  G(x,0)=G’(x,  l)=0. 
Thus  we  may  write  the  first  term  in  the  braces  in  (5.1.24) 


, 

>  j  d .  ■  • . 

y  •5','ifi  >  y}  ::  I  •t®sjosm£'Hioto 

■  '  :•  y  i;;.  j  ;•  'v1:,  ;/ 

■J  .a. -''O'. 

^•lv  ■  '  '  (  '  -!  i  ■  '  '  -r  at:  tcr-lrr,  (IS.  I,  r )  '  ■  rf 

■  .  ,  ■  .  ■ 

■  1  '>£  '  b i  OP"  O •  : BO  oW 

•  ( S--  .  >  .  “  '  I'D  .r  .'iCO  ..  X  '•  ‘  I C  .  ■  J. ■  ■;  l£| 


p.  p  ?,.i  o  •'  :j  ...  o/:  :  10  op- 


.  ■'  ,.j n  i:  v  p  ", '  '  o 
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as 

11m  G(x,h)  11m  G(x,h)  -C(x,0) 

h-0  h  “  l>-0  h 

(5.1.25) 


Hm  G(x,£ )  -G(x,  O)  dG(x,g) 

£-0  i  "  dC 


=  a  .  (5.1.26) 

By  a  similar  argument,  the  second  term  in  the  braces 
becomes 


-slim  G(x,  1-h)  _  pSG(x,  l) 

ph-0  h  ~  p  d£ 


(5.1.26) 


Thus,  we  may  write 

y(x)  -  p  .  a  st^-,01  k(?)d|. 

(5.1.27) 

Formula  (5. 1.27)  is  a  special  case  of  Green's 
formula.  In  passing,  we  remark  again  that  examination 
of  numerical  processes  may  often  point  toward  a  simple, 
intuitive  approach  to  theorems  of  analysis  (cf.  our 
discussion  of  the  fundamental  theorem  of  Calculus,  in 
Section  2.4). 
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The  preceeding  has  been,  at  best,  heuristic.. 

We  have  provided  no  proof  of  the  validity  of  our  limiting 
argument,  but  we  have  nevertheless  achieved  some  limited 
results.  In  the  following  section,  we  attempt  to  put  these 
results  on  a  sounder  footing. 

Section  5.2  A  Proof  of  Convergence 

In  solving  the  problem  (5.1.l),  we  replaced 
the  differential  equation  by  a  matrix  equation  to  obtain 
an  approximate  solution  £,  satisfying 

A  l  =  d .  (3.6.14) 

The  exact  solution  y_,  satisfies,  say 

A  ^  =  d  +  £.  (5.2.1) 

We  write 

z  =  z  =  n,  (5.2.2) 

and  hence,  substituting  into  (3.6,14),  we  obtain 

A  a  -  e .  (5.2.3) 

Clearly,  we  would  like  to  obtain  a  bound  on  iq 
and  show  that  -*■  0  as  h. 0,  To  do  so,  we  first  need  a 
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bound  on  e . 


In  deriving  (3.6.14),  we  used  the  approximation 


h2yi  *  62yi  =  yi+l-2yi+yl-l 


(5-2.4) 


Now,  by  Taylor’s  theorem. 


i+i  =  yi+hyi+lr  yi  +  fr  yi"’ +  sr  yiV(5ih  (5.2.5) 


h. 


4 

Yi  iv. 


-  2y±=  ~2y1> 


(5.2.6) 


and, 


,  ’  ,  h  ii 

yi-i  =  yi~hyi+  2T  yi 


-  yr  y^1  +  4Y  y^v  •  (5.2.7) 


4 

hH  iv 


Adding  these  three  together,  we  have 


e2yi  -  yi +  5i  (yiv(?i)+yiv(?2h 


(5.2.8) 


Assuming  that  y  is  analytic,  its  derivatives  are 
therefore  bounded,  and  we  may  find  some  number  7  satisfying 


5T  lyiV(lh)+yiV(|;2)  l<  y 


iv. 


(5.2.9) 


"  / 


o  no  6fiirocf 


.  d  ,  ::.v  .  v.  r  n  : 


,  -  ' 


c-'  ;  ..;T  v;,;  ;; 


q 


.  .  j  "■  ■■  ;■  ■ 
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uniformly  for  x  e[0,l],  or  for  i  e[0,n] .  Thus,  from  (5.2.8), 


62y1  =  hV  +  e 


(5-2.10) 


where,  for  each  element  of  £_,  we  have 


i  4 

£,  |<  7h  • 


(5.2.11) 


Now,  from  (5.2.3)* 


21  =  A  1e  , 


(5.2.12) 


and  therefore,  for  any  element,  r^,  of  tj_. 


n-1 


-1 


=  L  Ai^£r 

i=i 


(5.2.13) 


Thus , 


i=l 


<  yh 


4 


n-l 


i=l 


A 


-1 

i,  J0 1  * 


(5.2.14) 


Writing  the  expansion  of  A  in  terms  of  its 


eigenvalues  and  eigenvectors, 


-  't  i.  v  : ::  .  I  ?;  :  - 


,:i  .  -  -  lo.  .9  -i  c 


.  c'ir‘ 


. 

v-v ■•ocf  sjsnis-  n 
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A 


n-1 


I 

m  =1 


t 

m' 


(5.2.15) 


where 


=  X 


P  > 
m-m' 


and 


t 

EjR 


J 


(5.2.16) 


Therefore , 


°t  iM 

m=l  m 


(5.2.17) 


and  so 


Xmpm,ipm,^ ' 


(5.2.18) 


Accordingly,  from  (5.2.14), 


n-1 


< 7h  I 


m=l 


n-1 


i 


(5.2.19) 


hut,  it  is  well-known  that 


n-1 


<i  slh~l  <  \T~n  =  l/sfh, 


(5.2.20) 


I™/ 


...  :  ,6'.,  ?r{\ 


■  .  '  ,  1  .  i  t  4‘.;i  . 


/-■“  n 
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Therefore, 


(5-2.21) 


Suppose  that  the  X^  are  numbered  so  that  X^  Is  the  smallest 
in  absolute  value;  then 


<  _ L  =  *Y  h.3  .  (5.2.22) 

sTh  |\1|  sTh  |xj 

Thus,  we  see  that  unless  X^-+  0  faster  than  h.  , 

we  will  have  -*•  0,  as  h  -♦  0,  It  is  not  obvious  whether 

this  is  restrictive  or  not;  to  get  a  clearer  idea,  we  try 

an  example.  Let  us  suppose  that  in  (5.1.1),  q(x)  =  0.  Then, 

2 

the  matrix  A  becomes  the  well-known  matrix  D  ,  of  which  the 
smallest  eigenvalue  is  [29] 

=  4  sin2(2=r)  (5.2.23) 


(W  +  m 


nhx5 


}2 

•  •  •  J  , 


(5.2.24) 


j  j,  y 
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and  for  h  — ►  0 


2 


(5.2.25) 


Hence,  for  this  case 


(5.2.26) 


This  latter  result  encourages  us  to  believe 
that,  for  most  reasonable  differential  operators,  convergence 
should  be  obtained  with  little  difficulty. 

Henrici  [9]  proves  what  amounts  to  a  different, 
and  somewhat  weaker,  formulation  of  this  same  result,  in  a 
form  bearing  more  directly  on  the  discussion  of  the  preceeding 
section.  We  state  his  result,  modified  for  our  notation. 

If  g (x)  is  twice  continuously  differentiable  in 
[0,1],  andg(x)<  0,  then  there  exists  some  constant  K  such 
that 


a:1  -  G(ih, jh) k  h2K 
J 


(5.2.27) 


for  h  sufficiently  small.  Briefly,  this  implies  that  the 
correspondence  which  we  drew  between  equations  (5.1.21)  and 
(5.1.22)  is  valid. 


i  ■ 


.1  ,  ■  r.V  .  . 
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Section  5.3  An  Analytic  Expression  for  the  Green1 3 

Function 

In  practice,  as  already  indicated,  the  sytem 
(5.1.1)  subject  to  (5.1.2)  is  solved  by  replacing  it 
with 

A  £  =  d  =  h2k-b  (3.6.14) 

and  using  standard  numerical  techniques  to  solve  these 
linear  equations. 


In  this 

section. 

we  will 

exploit 

the  tri -diagonal 

feature 

of  A  to  obtain  an 

algebraic 

expression  for  A 

-1 

y 

and  hence,  to  G. 

Thence, 

we  will  obtain  an 

expression 

for  lim 

b  — ►  0  of 

G. 

Writing 

e  . 

1 

=  h2g1. 

(5.3.1) 

equation 

(5.1.12) 

becomes 

S 

V. 

1 — 1 
U) 

1 

CVJ 

)  -1 

0 

0 

•  •  • 

-1 

(2-e 

2)  -1 

0 

•  •  • 

0 

-1 

(2-e3) 

-1 

•  •  • 

-  A  = 

-1 

(2-e  , 

x  m-1 

)  -1 

(5.3. 

-1 

(2-e  ) 
x  my 

3 

'-vcn 


!  ■ 


fit 


;  '  xx  ,i  :  3;J< 
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where 


m  =  n-1. 


(5.3.3) 


We  define  AQ  =  1#  and  ^ (a^ ,a2, . . . ,a^)  to  be  the  deter¬ 
minant  of  an  i  by  l  matrix  of  the  form  of  (5.3.2)  with 
(2-ea  )  in  the  first  row,  (2-ea  )  in  the  second  row,  and 
so  on.  Then,  we  write  the  inverse  of  -A  in  terms  of  its 
co-factors.  For  m=2,  we  obviously  have  (since  G  =  h.A  1 

S 

A 


r  -  ~h 

-  -  a2u,2) 


Ax(2) 


A 


o 


A],  ( 1 ) 


\ 


(5.3.4) 


for  m  =  3,  we  have 

/ 


r  _  -h. 

-  A3(l,2,3) 


A2(2,3) 


A 


Ax(3) 


A^l) 


A, 


Ax(3)  a2(i,3)+a0  a1(i) 


A2(l,2) 


(5.3.5) 


The  elements  of  G  increase  rapidly  in  complexity 
as  m  increases,  but  by  taking  several  more  examples,  a 


b  :■■{.;  .  :  '  a  1 

,S *•/-.  ■ 


rf  i  /’ 


-i 
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pattern  emerges.  For  m=4,  we  find 


(1,2,34 


A3(2,3,4) 

A2(3,4) 


Ax(4) 


A2(3,4)  A1(4)  Aq 

a3(1;3,4)+a1(4)  a2(i^+aq  a1(i) 

A2(1;4)+Aq  A3(l,2;4)+A1(l)  A2(l,2) 

A1(l)  A2(l,2)  A3 (l,2,3) 

(5.3.6) 


given  by 


It  is  clear  from  these  examples  that  G^  ^  is 


°il  "  Am-i(i+1,i+2* • • °>m) 
h.  ( 1 , 2 ,  .  •  .  ,m) 


(5.3.7) 


and  G.^  by 


G. 


1  ,m 

h 


-Ai-l(i“1" i> • * • ^m_1) 

A^ (l,2, . . . ,m) 


(5.3.8) 


for  any  value  of  m.  In  a  similar  fashion,  we  find  the 
general  formula 


Am(l,2, . . .,m) 


h. 


G.  . 
J 


A 


>• 

■ 
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min£jn-i+l,  j) 

j^_j_  J  —  (2p-l)^^"',^,  *  *  *  *  J~  P  ■>  f  ~^"P  ■>  l^'P-t"l  )  •  •  M^l)) 


P=1 


(5.3.9) 


for  j  £  1, 
and 

oiJ  =  aJ±,  (5.3.10) 


for  j  >  i. 


We  now  require  an  expression  for  A^a^a^, .  .  „  <,a^) 
in  terms  of  the  £..  Consider  the  expansion  of  A„  in  terms 


i 


of  its  last  row 


Ai(al,a2'  *  *  *  ,CLll)  ~  eaJA£-l(ai’a2’°°*’CtJl-l') 


A£- 2^al,a2’  *  * ‘ ,a£-2^  * 


(5.3.11) 


Repeating  this  procedure,  we  attempt  to  write  A^  as  a 


polynomial  of  the  form 


A^(ai,a2, . . .,a^)  -  l+l  +  \  a^p  £a 


l 

l 

P=1 


_  _  a^.Pi>p2e“p  e% 

p2=l  Pl=l  d  P1  P2 


P 


P, 


:  J  s  ,  ■  '3  .. D  •  3  c ■  3'  q  ?  .c  £  f ::  ! j: ' ;  £  ;•  q  -j  ft 

' 
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l  p3  P2 


+ 


L  L  L  ^Pi  >  Pn '  Pq  an  an 
=  1  P  =1  P  =1  X  ^  D  ^1  ^ 


£~  £rv  +  •  •  •  > 

\JL 


(5.3.12) 


P3=i  p2=i  pr 


2  p3 


where  the  a*s  are  still  to  be  determined.  By  comparing 
(5.3.12)  with,  the  expansion  for  the  determinant  we 
find 


i 

^  f  ^2  *  *  *  *  *  H  )  =  ~  (^il~P)P  £  ££ 

P=1  P 


i  P2 


I  I  (^+1-P2)(P2-p1)pl  £an  £ 


a 


p2=1  P^l 


P]_  ^2 


^  P3  p2 


'  I  I  I  (i+l-p3)(p3-p2)(p2-p1)p1ea>  Ba  E+a; 


P3=l  p2=l  P]_=l 


PX  P2  P3 


(5-3.13 


We  now  -  in  view  of  later  applications  -  consider 

the  special  case  of  ^=m  of  (5.3.13).  We  have 

m 

hAm(l j 2, . .  .,m)  =  (m+l)h  -  V  [(m+l)h-ph.]phg(ph.)h 

P=1 

m  p2 

+  Z  Z  ^+1)^'~P2h^  tp2h“pih] pihs(piln)g(p2h)hh 

p2=l  P^l 


+ . . . 


(5.3.14) 


'  .  j.  .  ;r  :scf ...-  rj:  S  £  ,  £  V) 


j. 


,  ''  1  '.  *  ■  • 


rt 


a  i,:  !i;  ViJ ' '  , 


1  .2 


Now,  setting 


u.  p.h  =  p  /(m+l) 


(5 


we  have 


1-h 


hAm(l,2, . . . ,m)  =  1-  ^  [l-u]ug(u)h 


u=h 


1-h  u2 


+ 


V  [ 1 -Up ] [u2-u^ ]u1g(u1)g(u2)hh+, 


u2=h  ufh 


(5 


h  — ►  0  u=0 


r  r 

1  -  j  (l-u)ug(u)du 


n  —  n  u2 

+  J  J  ^(l-u2)  (u2-u-L)u1g(u1)g(u2)du1du2+. 
u2=0  u-L=0  ^ 


We  must  now  consider  another  special  case, 
first  column  of  the  inverse  matrix.  We  have,  from 

A. 

TT  Gil  =  i+2,  . .  .  ,m) , 


3.15) 


3.16) 


3.17) 


the 

(5.3.7) 


and  hence,  using  (5.3.13) 


«JXI 
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m-i 

-  A.  '  (m+l)h-ih- V  [  (m+l)h-ih-ph]phg(ph+ih)h 

P=1 


g(p2h+ih)h.h- .  .  . , 


h 


*  1-x 

0 


-x)ug(u+x)du 


(5.3.18) 


1-X  n  u. 


+ 


u2-0  u^=0 


(l-uQ-x) (u0-uu )un g(un  +x)g(u0+x)du1  du, 


which,  we  write  as  $(x). 


Clearly  then,  from  (5-3.17) , 


(5.3.19) 


hAm(l,2,...,m)-  4>(0), 


(5.3.20) 


and  so  we  write 


G.  hA 

1,1  _  _ m - 1  <3Hx 


h  hAm  0  ( 0 


=  £}(x) 


(5.3.21) 


1 
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Evidently, 

n(o)  =  |[2|  =  1,  (5.3.22) 

and  also,  from  (5-3.19) >  since,  for  x  =  1,  the  Integrals 
are  over  a  range  of  zero  length. 


Q(l)  =--  0. 


(5.3.23) 


It  is  fairly  clear  that,  by  a  process  similar 
to  that  just  described,  we  will  also  have 


(5.3.24) 


From  (5.1.24)  and  (5.1.27),  we  see  that  fi(x) 
from  (5.3.21)  Is 

fl(.x)  =  ^§^1  ’  (5.3.25) 

and  that  ^(l-x),  from  (5-3-24)  Is 

n(l-x)  =  -  (5.3.26) 

Therefore,  we  see  that 

a^(x)  =  n(x),  (5.3.27) 
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and 


Oi^x)  =  fi(l-x). 


(5.3.28) 


are  two  independent  solutions  of 


cd"  +  goo  =  0 


(5.3.29) 


subject  to 


oXj_(0)  =  av,(l)  =  1 


(5.3.30) 


and 

0^(1)  =  o>2(0)  =  0.  (5.3.31) 


Hence ,  according  to  a  standard  result  [28]  we 
can  now  write 


f~ I  to2(^)^1(l)  >  x<  i 

G(x^)  =  j  W  (5.3.32) 


where  W  is  Wronskian 

t  t 

W  =  CD^OD^  -  ^2^1 


(5.3.33) 


and  is  a  constant.  We  evaluate  W  at  x  =  1,  and  hence 
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(5.3.34) 


but  $  (l)  Is  easily  seen  to  be  -1,  and  so 


W  =  =  l/$(o) . 


(5.3.35) 


We  thus  rewrite  (5.3.32) 


(5.3.36) 


where  $(x)  is  given  by  (5.3.19). 

The  line  of  argument  throughout  this  section  has 
not  been  rigorous;  in  particular.,  we  have  provided  no 
guarantee  that  the  limits  we  have  written  actually  exist. 

We  gain  some  confidence  (but  not  proof)  from  the  theorem  of 
the  preceeding  section;  we  will  attempt  further  to  make 
the  argument  seem  plausible  by  means  of  the  examples  of 
the  next  section. 

Section  5.4  Special  Cases  of  the  General  Formula 

We  shall  verify  only  that  the  cd!  s  defined  by 


(5.3.27)  and  (5.3.28)  actually  satisfy  (5.3.29)  to  (5*3.31) 


% 
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In  some  particular  Instance since  the  work  beyond  that 
point  is  all  well-known.  It  would  be  more  satisfying  to 

complete  the  anlytical  investigation  of  the  difficult 
limit  processes  involved,  but  this  does  not  appear  tractable 
at  this  stage. 


First,  we  try  g(x)=0.  Then,  we  have 


(5-4.1) 


and  evaluating  (5.3.19)  with  g  =  0,  all  the  integral 


terms  vanish  and  we  are  left  with 


<t>(x)  =  1-x, 


(5.4.2) 


and  hence 


4>(0)  =  1. 


(5.4.3) 


We  test  =  1-x  in  the  differential  equations: 


(5.4.4) 


and 


(5.4.5) 


and  so  (5.3.29)  is  satisfied.  Furthermore. 
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IX,  (0)  1,  and  av,  ( 1 )  0,  (5.4.6) 

and  so  (5*3. 30)  and  (5-3.31)  are  satisfied. 


Also, 

0)o(x)  =  ft(l-x)  =  $(l-x)/$(o) 

=  l-(l-x)  =  X  (5-4.7) 

and,  as  before,  the  required  conditions  are  satisfied. 

We  now  try  the  slightly  more  difficult  function, 
g(x)  =  1.  Thus,  (5.3.29)  becomes 


oy'+cD  -  0 


(5.4.8) 


Now,  from  (h.3.19);  with  g  =  1,  we  have 


p  1  x 

<E>(x)  =  1-x  -  /  (l-u-x)u  du 


1-X 


+ 


(l_Up-x)  (^-u-^u-^du^dug-  ...  .  (5.4.9) 


Up=0  u^=0 


The  first  integral  term  is 


1-x 


1-x 


h  = 


u=0 


(l-u-x)u  du  =  (l-x)  /  udu  - 


1-x 


2 

u  du 


u=0 


u-0 
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1-X 

(1-x) 

u 

L  2  - 

0 

(1-x)3  

Liz. 

u 


3‘ 


1-x 


3 


Jo 


3 


3 


(5.4.10) 


The  second  integral  term  is 


1-x 


u. 


In  =  /  (l-uQ-x)  /  (uQ-un )un dun duQ 
^u2=0  d  J  UpO  d  1  1  1  ^ 


u20 


3 

-Q 

(l-u2-x)  du,„=  ( 1-x) 


1-x 


3 


u 


n  1-X. 

_2  du0-/  u^du0 

(C  2  J  2d 

u2=0 


u=0 


(l-x)-:  (1-x)5  _  (l-x)5  _  ( 1  — x) 5 

25  30-  120  55 


(5.4J1) 


It  can  be  shown  that 


i  _^\3  f1  _v\5 


®(x)  =  (l-X)-^-  + 


3 


55 


(5.4.12) 


which,  we  recognize  as  the  Taylor  series  for  sin(l-x) 
hence 


$(x)  =  sin(l-x) ; 
$(o)  =  sin(.l); 
i>(l-x)  =  sin(x)  . 


(5.4.13) 


J 


’l 
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Thus , 


to-^x)  =  sin(l-x)/sin(l) ,  (5.4.14) 


and 


oi->(x)  =  sin(x)/sin(l) . 


(5.4.15) 


Now, 


cUj_(x)  =  -  cos (l-x)/sin(l) ,  (5.4.16) 


and. 


(i^(x)  =  -  sin(l~x)/sln(l) ,  (5.4.17) 


and  so. 


“l  +  ^  =  0 


(5.4.18) 


and  (5.3.29)  is  satisfied.  Clearly 

0^(0)  =  sin(l)/sin(l)  =  1,  (5.4.19) 

and 

515(1)  -  sin(0)/sin(l)  =  0,  (5.4.20) 

satisfying  (5.3.30)  and  (5.3.31). 

The  details  for  u0(x)  are  equally  trivial.  We 
have  thus  obtained  concrete  verification  of  the  formulas 
of  the  previous  section,  and  can  feel  fairly  confident 


■ 


analysis  performed  there  is  correct,  and  will 
probably  work  for  more  complicated  functions,  g(x), 
although  the  integrations  may  become  very  difficult  to 
perform.  It  may  be  that  a  possible  approach  is  to 
approximate  g(x)  by  a  polynomial,  and  then  appeal  to  the 
Weierstrass  approximation  theorem,  but  even  polynomials 
render  the  integrations  in  (5.3.19)  very  difficult. 

It  is  unlikely  that  the  formulas  which,  we 
have  found  are  in  any  way  practical  devices  for  finding 
analytic  solutions  to  equations.  They  may,  however,  have 
some  value  in  finding  numerical  solutions. 
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CHAPTER  VI 

CONCLUDING  REMARKS 

In  the  preceeding  chapters,  we  have  reviewed 
a  narrow,  but  important,  spectrum  of  numerical  techniques 
for  solving  differential  equations,  and  have  discussed  in 
some  detail  certain  difficulties  -  such,  as  instability  - 
which  hamper  us  in  our  attempts  to  obtain  greater  precision 
in  the  numerical  results.  Partial  differential  equations 
have  been  treated  very  briefly,  but  at  sufficient  length 
to  show  that  each  of  the  three  main  types  requires  a 
different  numerical  approach  -  one  which  is  suggested  by 
the  analytic  character  of  the  differential  system  itself. 

In  one  chapter,  wre  have  examined  the  analytical  behaviour 
of  the  finite -difference  analogue  of  an  ordinary  dif¬ 
ferential  equation  subject  to  boundary  conditions,  as  the 
tabular  interval,  h,  approaches  zero.  The  analysis  is 
incomplete  at  present,  and  may  be  regarded  as  being  of 
mostly  heuristic  value. 

The  numerical  study  of  differential  equations  may 
seem  to  the  mathematician  nothing  but  a  hodge-podge  of 
special  devices.  It  may,  however,  be  said  with  some 
confidence  that  progress  has  now  brought  us  to  a  point  at 
which  very  few  ordinary  differential  equations  present  any 
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serious  difficulties  to  the  numerical  analyst.  Probably, 
little  more  than  this  should  be  expected  at  this  time, 
because  the  intensive  study  of  the  numerical  solution  of 
ordinary  differential  equations  began  only  recently,  and 
because  it  has  usually  been  driven  by  the  day-to-day 
necessity  of  solving  this  or  that  particular  problem. 

It  need  hardly  be  emphasized  that  this  situation  Is 
unsatisfactory. 

Von  Neumann  once  remarked  that  large  digital 
computers  could  be  justified  simply  on  the  basis  that 
numerical  studies  of  non-linear  differential  systems  might 
lend  some  insight  into  those  topics  which  ought  to  be 
studied  by  an  analyst;  or,  as  Hamming  has  said  more 
recently,  "The  purpose  of  computing  is  insight,  not  numbers 
Von  Neumann* s  hope  has  not  yet  been  fulfilled,  but  his 
remark  does  explain,  and  to  a  large  extent  justify,  the 
energy  expended  by  numerical  workers  on  particular  problems 
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